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Abstract 

Explicit  expressions  for  the  stiffness,  geometric  stiffness,  mass,  and  aerodynamic  force  matrices 
are  derived  for  the  flutter  analysis  of  simply  supported  composite  wing  panels.  The  formulation  is  based 
on  Ritz  analysis  using  simple  polynomials  and  Piston  Theory  aerodynamics.  The  use  of  simple 
polynomials  eliminates  the  need  for  numerical  quadrature  and  decreases  computation  time.  Analytic 
sensitivities  of  the  aeroelastic  system  matrices  and  critical  dynamic  pressures  are  obtained  with  respect  to 
layer  thickness,  fiber  direction,  and  panel  shape.  The  method  is  integrated  with  wing  box  analysis  based 
on  either  the  equivalent  plate  approach  or  finite  element  method,  making  it  possible  to  obtain  sensitivities 
of  the  stability  boimdary  with  respect  to  wing  planform  shape  or  locations  of  ribs  and  spars.  The  analytic 
sensitivities  are  used  to  construct  approximations  of  the  aeroelastic  stability  boundary  for  integrated  wing  / 
panel  design  synthesis. 


ABBREVIATED  BIBLIOGRAPHY 


1.  Dowell,  E.H.,  Aeroelasticity  of  Plates  and  Shells,  Noordhoff  International  Publishers,  Leyden,  1975. 

2.  Sander,  G.,  Bon,  C.,  and  Geradin,  M.,  “Finite  Element  Analysis  of  Supersonic  Panel  Flutter”, 
Intemation  Journal  for  Numerical  Methods  in  Engineering,  7,  pp.  379-394,  1973. 

3.  Bismarck-Nasr,  M.,N.,  “Finite  Element  Analysis  of  Aeroelasticity  of  Plates  and  Shells”,  ASME 
Applied  Mechanics  Reviews,  Vol.  45,  No.  12,  pp.  461-482,  December  1992. 

4.  Hajela,  P.,  and  Glowasky,  R.,  “Application  of  Piezoelectric  Elements  in  Supersonic  Panel  Flutter 
Suppression”,  AIAA  Paper  91-3191,  AIAA  Aircraft  Design  Systems  and  Operations  Meeting,  Baltimore, 
MD,  September  1991. 

5.  Scott,  R.C.,  and  Weishaar,  T.A.,  “Controlling  Panel  Flutter  Using  Adaptive  Materials”,  AIAA  Paper 
91-1007,  April  1991. 

6.  Dowell,  E.H.,  and  Ilgamov,  M.,  Studies  in  Nonlinear  Aeroelasticity,  Springer- Verlag,  New  York,  1988. 

7.  Frampton,  K.D.,  Clark,  R.L.,  and  Dowell,  E.H.,  “State  Space  Modeling  for  Aeroelastic  Panels  Subject 
to  Linearized  Potential  Flow  Aerodynamic  Loading”,  AIAA  Paper  95-1294,  Proceedings  of  the  36th 
AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics,  and  Materials  Conference,  New 
Orleans,  LA,  April  10-13,  1995,  Vol.  2,  pp.  1183-1189. 

8.  Kariappa,  Somashekar,  B.R.,  and  Shah,  C.G.,  “Discrete  Element  Approach  to  Flutter  of  Skew  Panels 
with  IN-Plane  Forces  Under  Yawed  Supersonic  Flow”,  AIAA  Journal,  Vol.  8,  No.  11,  November  1970, 
pp.  2017-2022. 

9.  Chowdary,  T.V.R.,  Sinha,  P.K.,  and  Parthan,  S.,  “Finite  Element  Flutter  Analysis  of  Composite  Skew 
Panels”,  Computers  &  Structures,  Vol.  58,  No.  3,  pp.  613-620,  1996. 

10.  Voss,  H.  M.,  and  Dowell,  E.  H.,  “Effect  of  Aerodynamic  Damping  on  Flutter  of  Thin  Panels”,  AIAA 
Journal,  Vol.  2,  No.  1,  January  1964,  ppl  19-120. 

11.  Lottatti,  I.,  “The  Role  of  Damping  on  Supersonic  Panel  Flutter”,  AIAA  Journal,  Vol.  23,  No.  10, 
October  1985,  pp.  1640-1642. 

12.  Giles,  G.L.,  “Equivalent  Plate  Analysis  of  Aircraft  Wing  Box  Structures  with  General  Planform 
Geometry”,  Journal  of  Aircraft,  Vol.  23,  No.  11,  November  1986,  pp.  859-864. 

13.  Giles,  G.L.,  “Further  Generalization  of  an  Equivalent  Plate  Representation  for  Aircraft  Structrrral 
Analysis”,  Journal  of  Aircraft,  Vol.  26,  No.  1,  Janrtary  1989,  pp.  67-74. 

14.  Livne,  E.,  “Integrated  Multidisciplinary  Optimization  of  Actively  Controlled  Fiber  Composite  Wings”, 
Ph.D.  dissertation.  Mechanical,  Aerospace  and  Nuclear  Engineering  Department,  University  of 
California,  Los  Angeles,  1990. 

15.  Livne,  E.,  “Analytical  Sensitivities  for  Shape  Optimization  in  Equivalent  Plate  Structural  Wing 
Models”,  Journal  of  Aircraft,  Vol.  31,  No.  4,  pp.  961-969,  July-August  1994. 

16.  Livne,  E.,  and  Li,  W-L.,  “Aeroservoelastic  Aspects  of  Wing  /  Control  Surface  Planform  Shape 
Optimization”,  AIAA  Journal,  Vol.  33,  No.  2,  February  1995,  pp.  302-311. 


17.  Livne,  E.,  and  Milosavljevic,  R.,  “Analytic  Sensitivity  and  Approximation  of  Skin  Buckling 
Constraints  in  Wing  Shape  Synthesis”,  Journal  of  Aircraft,  Vol.  32,  No.  5,  September-October  1995,  pp. 
1102-1113. 

18.  Li,  W-L.,  and  Livne,  E.,  “Analytic  Sensitivities  and  Approximations  in  Supersonic  and  Subsonic 
Wing  /  Control  Surface  Unsteady  Aerodynamics”,  AIAA  Paper  95-1219,  AIAA/ASME/ASCE/AHS/ASC 
36th  Structures,  Structural  Dynamics  and  Materials  Conference,  New  Orleans,  April  1995. 

19.  Haftka,  R.T.,  and  Gurdal,  Z.,  “Elements  of  Structural  Optimization”,  Third  Edition,  Kluwer  Academic 
Publishers,  Dordrecht,  Boston,  London  1991.  (Chapter  6) 

20.  Pedersen,  P.,  and  Seyranian,  A.,  “Sensitivity  Analysis  for  Problems  of  Dynamic  Stability”, 
International  Joimial  of  Solids  and  Structures,  Vol.  19,  No.  4,  pp.  315-335,  1983. 

21.  Seyranian,  A.,  “Sensitivity  Analysis  of  Multiple  Eigenvalues”,  Mechanical  Structures  and  Machines, 
Vol.  21,  No.  2,  pp.  261-284,  1993. 

22.  Golub,  G.,  and  Van  Loan,  C.  F.,  Matrix  Computations,  2nd  edition,  Johns  Hopkins  University  Press, 
1989. 

23.  Lecina,  G.,  and  Petiau,  C.,  “Advances  in  Optimal  Design  with  Composite  Materials”,  in  Computer 
Aided  Optimal  Design;  Structuyral  and  Mechanical  Systems,  C.A.  Mota  Soares,  editor.  Springer  Verlag, 
Berlin,  1987,  pp.  943-953. 

24.  Ausman,  J.D.,  Hangen,  J.A.,  and  Acevedo,  D.A.,  “Application  of  a  Local  Panel  Buckling  Constraint 
within  Automated  Multidisciplinary  Structural  Analysis  and  Design”,  AIAA  92-1116,  1992  Aerospace 
Design  Conference,  February  3-6,  1992,  Irvine,  CA. 

25.  Bruhn,  E.F.,  Analysis  and  Design  of  Flight  Vehicle  Structures,  Tri  State  Offset  Company,  1973, 
Chapter  C5. 

26.  Whitney,  J.  M.,  Structural  Analysis  of  Laminated  Anisotropic  Plates,  Technomic  Publishing 
Company,  Lancaster,  PA,  1987. 

27.  Tsai,  S.W.,  and  Hahn,  H.T.,  Introduction  to  Composite  Materials,  Technomic  Publishing  Company, 
Lancaster,  PA,  1980. 


FLUTTER  ANALYSIS  AND  ANALYTIC  SENSITIVITIES 
FOR  TRAPEZOIDAL  PANELS 


by 

David  Mineau 
2iid  Lieutenant,  USAF 


Abstract 

Explicit  expressions  for  the  stiffness,  geometric  stiflhess,  mass,  and  aerodynamic  force  matrices 
are  derived  for  the  flutter  analysis  of  simply  supported  composite  wing  panels.  The  formulation  is  based 
on  Ritz  analysis  using  simple  polynomials  and  Piston  Theory  aerodynamics.  The  use  of  simple 
polynomials  eliminates  the  need  for  numerical  quadrature  and  decreases  computation  time.  Analytic 
sensitivities  of  the  aeroelastic  system  matrices  and  critical  dynamic  pressures  are  obtained  with  respect  to 
layer  thickness,  fiber  direction,  and  panel  shape.  The  method  is  integrated  with  wing  box  analysis  based 
on  either  the  equivalent  plate  approach  or  finite  element  method,  making  it  possible  to  obtain  sensitivities 
of  the  stability  boundary  with  respect  to  wing  planform  shape  or  locations  of  ribs  and  spars.  The  analytic 
sensitivities  are  used  to  construct  approximations  of  the  aeroelastic  stability  boundary  for  integrated  wing  / 
panel  design  synthesis. 


ABBREVIATED  BIBLIOGRAPHY 


1.  Dowell,  E.H.,  Aeroelasticity  of  Plates  and  Shells,  Noordhoff  International  Publishers,  Leyden,  1975. 

2.  Sander,  G.,  Bon,  C.,  and  Geradin,  M.,  “Finite  Element  Analysis  of  Supersonic  Panel  Flutter”, 
Intemation  Journal  for  Numerical  Methods  in  Engineering,  7,  pp.  379-394,  1973. 

3.  Bismarck-Nasr,  M.,N.,  “Finite  Element  Analysis  of  Aeroelasticity  of  Plates  and  Shells”,  ASME 
Applied  Mechanics  Reviews,  Vol.  45,  No.  12,  pp.  461-482,  December  1992. 

4.  Hajela,  P.,  and  Glowasky,  R.,  “Application  of  Piezoelectric  Elements  in  Supersonic  Panel  Flutter 
Suppression”,  AIAA  Paper  91-3191,  AIAA  Aircraft  Design  Systems  and  Operations  Meeting,  Baltimore, 
MD,  September  1991. 

5.  Scott,  R.C.,  and  Weishaar,  T.A.,  “Controlling  Panel  Flutter  Using  Adaptive  Materials”,  AIAA  Paper 
91-1007,  April  1991. 

6.  Dowell,  E.H.,  and  Ilgamov,  M.,  Studies  in  Nonlinear  Aeroelasticity,  Springer- Verlag,  New  York,  1988. 

7.  Frampton,  K.D.,  Clark,  R.L.,  and  Dowell,  E.H.,  “State  Space  Modeling  for  Aeroelastic  Panels  Subject 
to  Linearized  Potential  Flow  Aerodynamic  Loading”,  AIAA  Paper  95-1294,  Proceedings  of  the  36th 
AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics,  and  Materials  Conference,  New 
Orleans,  LA,  April  10-13,  1995,  Vol.  2,  pp.  1183-1189. 

8.  Kariappa,  Somashekar,  B.R.,  and  Shah,  C.G.,  “Discrete  Element  Approach  to  Flutter  of  Skew  Panels 
with  IN-Plane  Forces  Under  Yawed  Supersonic  Flow”,  AIAA  Journal,  Vol.  8,  No.  11,  November  1970, 
pp.  2017-2022, 

9.  Chowdary,  T.V.R.,  Sinha,  P.K.,  and  Parthan,  S.,  “Finite  Element  Flutter  Analysis  of  Composite  Skew 
Panels”,  Computers  &  Structures,  Vol.  58,  No.  3,  pp.  613-620, 1996. 

10.  Voss,  H.  M.,  and  Dowell,  E.  H.,  “Effect  of  Aerodynamic  Damping  on  Flutter  of  Thin  Panels”,  AIAA 
Journal,  Vol.  2,  No.  1,  January  1964,  ppl  19-120. 

11.  Lottatti,  L,  “The  Role  of  Damping  on  Supersonic  Panel  Flutter”,  AIAA  Journal,  Vol.  23,  No.  10, 
October  1985,  pp.  1640-1642. 

12.  Giles,  G.L.,  “Equivalent  Plate  Analysis  of  Aircraft  Wing  Box  Structures  with  General  Planform 
Geometry”,  Journal  of  Aircraft,  Vol.  23,  No.  11,  November  1986,  pp.  859-864. 

13.  Giles,  G.L.,  “Further  Generalization  of  an  Equivalent  Plate  Representation  for  Aircraft  Structural 
Analysis”,  Journal  of  Aircraft,  Vol.  26,  No.  1,  January  1989,  pp.  67-74. 

14.  Livne,  E.,  “Integrated  Multidisciplinary  Optimization  of  Actively  Controlled  Fiber  Composite  Wings”, 
Ph.D.  dissertation.  Mechanical,  Aerospace  and  Nuclear  Engineering  Department,  University  of 
California,  Los  Angeles,  1990. 

15.  Livne,  E.,  “Analytical  Sensitivities  for  Shape  Optimization  in  Equivalent  Plate  Structural  Wing 
Models”,  Journal  of  Aircraft,  Vol.  31,  No.  4,  pp.  961-969,  July-August  1994. 

16.  Livne,  E.,  and  Li,  W-L.,  “Aeroservoelastic  Aspects  of  Wing  /  Control  Surface  Planform  Shape 
Optimization”,  AIAA  Journal,  Vol.  33,  No.  2,  February  1995,  pp.  302-311. 


17.  Livne,  E.,  and  Milosavljevic,  R.,  “Analytic  Sensitivity  and  Approximation  of  Skin  Buckling 
Constraints  in  Wing  Shape  Synthesis”,  Journal  of  Aircraft,  Vol.  32,  No.  5,  September-October  1995,  pp. 
1102-1113. 

18.  Li,  W-L.,  and  Livne,  E.,  “Analytic  Sensitivities  and  Approximations  in  Supersonic  and  Subsonic 
Wing  /  Control  Surface  Unsteady  Aerodynamics”,  AIAA  Paper  95-1219,  AIAA/ASME/ASCE/AHS/ASC 
36th  Structures,  Structural  Dynamics  and  Materials  Conference,  New  Orleans,  April  1995. 

19.  Haftka,  R.T.,  and  Gurdal,  Z.,  “Elements  of  Structural  Optimization”,  Third  Edition,  Kluwer  Academic 
Publishers,  Dordrecht,  Boston,  London  1991.  (Chapter  6) 

20.  Pedersen,  P.,  and  Seyranian,  A.,  “Sensitivity  Analysis  for  Problems  of  Dynamic  Stability”, 
International  Journal  of  Solids  and  Structures,  Vol.  19,  No.  4,  pp,  315-335,  1983. 

21.  Seyranian,  A.,  “Sensitivity  Analysis  of  Multiple  Eigenvalues”,  Mechanical  Structures  and  Machines, 
Vol.  21,  No.  2,  pp.  261-284,  1993. 

22.  Golub,  G.,  and  Van  Loan,  C.  F.,  Matrix  Computations,  2nd  edition,  Johns  Hopkins  University  Press, 
1989. 

23.  Lecina,  G.,  and  Petiau,  C.,  “Advances  in  Optimal  Design  with  Composite  Materials”,  in  Computer 
Aided  Optimal  Design:  Structuyral  and  Mechanical  Systems,  C.  A.  Mota  Soares,  editor.  Springer  Verlag, 
Berlin,  1987,  pp.  943-953. 

24.  Ausman,  J.D,,  Hangen,  J,  A.,  and  Acevedo,  D.  A.,  “Application  of  a  Local  Panel  Buckling  Constraint 
within  Automated  Multidisciplinaiy  Structural  Analysis  and  Design”,  AIAA  92-1116,  1992  Aerospace 
Design  Conference,  February  3-6,  1992,  Irvine,  CA. 

25.  Bruhn,  E.F.,  Analysis  and  Design  of  Flight  Vehicle  Structures,  Tri  State  Offset  Company,  1973, 
Chapter  C5. 

26.  Whitney,  J.  M.,  Structural  Analysis  of  Laminated  Anisotropic  Plates,  Technomic  Publishing 
Company,  Lancaster,  PA,  1987. 

27.  Tsai,  S.W.,  and  Hahn,  H.T.,  Introduction  to  Composite  Materials,  Technomic  Publishing  Company, 
Lancaster,  PA,  1980. 


FLUTTER  ANALYSIS  AND  ANALYTIC  SENSITIVITIES 
FOR  TRAPEZOIDAL  PANELS 


by 

David  Mineau 


A  thesis  submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

Master  of  Science 
in  Aeronautics  and  Astronautics 

University  of  Washington 

1996 


Approved  by _  el  _ 

(Chairperson  of  Supervisory  Committee) 

Program  Authorized 

to  Offer  Degree _ 


Date 


Master^s  Thesis 

In  presenting  this  thesis  in  partial  fulfillment  of  the  requirements  for  a  Master’s  degree  at  the  University 
of  Washington,  I  agree  that  the  Library  shall  make  its  copies  freely  available  for  inspection.  I  further 
agree  that  extensive  copying  of  this  thesis  is  allowable  only  for  scholarly  purposes,  consistent  with  “fair 
use”  as  prescribed  in  the  U.S.  Copyright  Law.  Any  other  reproduction  for  any  purposes  or  by  any  means 
shall  not  be  allowed  without  my  written  permission. 

Signature _ 

Date _ 


TABLE  OF  CONTENTS 


Page 

List  of  Figures .  ii 

List  of  Tables .  iii 

Nomenclature  .  iv 

Chapter  1:  Introduction .  1 

Chapter  2:  Panel  Flutter  Equations  of  Motion .  3 

Chapter  3:  Polynomial  Modeling .  16 

Chapter  4:  Mass  Matrix . 21 

Chapter  5:  Stiffness  Matrix .  23 

Chapter  6:  Geometric  Stiffness  Matrix .  28 

Chapter?:  Aerodynamic  Force  Matrices .  36 

Chapter  8:  Analytic  Sensitivities .  42 

Chapter  9:  Verification  of  Results .  56 

Chapter  10:  Conclusion . 74 

Bibliography . 76 


LISTOFHGURES 


Number  Page 

2.1  Stability  Analysis  in  the  Presence  of  Damping  .  14 

2.2  Stability  Analysis  in  the  Absence  of  Damping  .  15 

3.1  Wing  Planform  .  17 

3.2  Panel  Layout  .  18 

7. 1  Panel  in  Angled  Flow  .  36 

9.1  Convergence  of  Solution  for  Critical  Dynamic  Pressure  with  no  In-Plane  Loads  .  57 

9.2  Panel  Geometry  for  Ref.  5  .  57 

9.3  Panel  Geometry  for  Ref.  21  .  59 

9.4  Panel  Geometry  for  Ref.  23  .  61 

9.5  Acrit  versus  Fiber  Angle  for  Different  Skew  Angles  .  64 

9.6  iJcrit  versus  Fiber  Angle  for  Different  Skew  Angles  .  65 

9.7  Fiber  Angle  Sensitivity  vs.  Step  Size  . . .  66 

9.8  Acrit  vs.  Fiber  Angle  for  \|r=15°  . .  67 

9.9  Anutt  vs.  Fiber  Angle  for  \jr=15°  .  68 

9. 10  Acrit  vs.  Fiber  Angle  for  xjr=45°  .  69 

9.11  Anutt  vs.  Fiber  Angle  for  tj/=45°  . , .  70 

9.12  Panel  Geometry  For  Xar  Sensitivity  .  71 

9.13  Anut  vs.  Xar  .  72 

9.14  Anutt  vs.  In-plane  Load  Factor  . . .  73 


ii 


LIST  OF  TABLES 


Number  Page 

3. 1  Constants  Ui  and  their  corresponding  power  of  y  . .  17 

3.2  Constants  Vj  and  their  corresponding  power  of  x  and  y  .  18 

5.1  Coefficients  and  powers  of  polynomial  terms  in  the  [F3]  Matrix  .  26 

6.1  Coefficients  and  Powers  of  x  and  y  of  the  matrix  Alpha  .  31 

6.2  Coefficients  and  Powers  of  Terms  of  the  [F2]  Matrix  .  34 

9.1  Acrit  for  Rectangular  Isotropic  Panels  .  58 

9.2  Acrit  for  Skewed  Panels  with  (j)  =  0"^  . 59 

9.3  Acrit  for  Skewed  Panels  with  (j)  =  15"^  .  60 

9.4  Anutter  and  Acrit  for  Skewed  Panels  .  60 

9.5  Critical  Dynamic  Pressures  With  Tyy  =  Txy  =  0  .  61 

9.6  Critical  Dynamic  Pressures  With  r^x  =  Tyy  =  0  .  61 

9.7  Acrit  for  Various  Skew  Angles  .  62 

9.8  Acrit  for  Skewed  Composite  Panel  .  62 

9.9  Qcrit  for  Skewed  Composite  Panels  . 63 

9.10  Aflutter  and  Qnutter  for  Skewed  Composite  Panels  .  63 

9.11  Fiber  Angle  Sensitivities  at  0=45"^  . . . .  66 


iii 


NOMENCLATURE 


[  A]  in-plane  local  stiffness  matrix  for  the  plate  (3  x  3) 

\Adamp\  aerodynamic  damping  matrix 

\-^stiff\  aerodynamic  stiffness  matrix 

cj  total  system  damping  matrix 

[  D]  out  of  plane  local  stifftiess  matrix  for  the  plate  (3  x  3) 

/,  (t,  y)  admissible  functions 

[  ^1  ]  >  [  ^2  ]  ’  [  ^3  ]  matrices  containing  admissible  functions  and  their  derivatives 

T)  weight  function  for  admissible  functions  ensuring  zero  displacement  on  panel  boundary 


Kx,y) 

Hix,y) 

K 

I 5^2 

[K],[K^] 

P] 

m,n 

mi,ni 

mh,n^ 

V  V  M 


coefficients  in  the  polynomial  expression  for  [F3  ] 
total  thickness  of  panel 
depth  of  wing 

coefficients  of  wing  depth  polynomial 

indices  of  terms  in  layer  thickness  polynomials 

integral  of  a  simple  polynomial  term  over  trapezoidal  panel  area 

panel  stiffness  and  geometric  stiffness  matrix,  respectively 

total  system  stiffness  matrix 

powers  of  x  and  y  in  a  polynomial  term 

powers  of  x  and  y  in  polynomial  series  for  thickness  of  layer  i 

powers  of  x  and  y  in  Ritz  functions 

powers  of  x  and  y  in  the  wing  depth  series 

powers  of  x  and  y  in  polynomial  terms  of  [F3  ] 


mf  2f  ,nf  2  f  powers  of  x  and  y  in  elements  of  [ F2  ] 


rrtqw.pw  ,nqw,pw  powers  of  x  and  y  for  elements  of  [a] 


panel  mass  matrix  and  total  system  mass  matrix 


IV 


N. 

i 

N, 

[N] 


free  stream  Mach  number 

number  of  terms  in  the  thickness  polynomial  for  layer  i 

number  of  layers 

2x2  matrix  of  in-plane  loads 

in-plane  loads  per  unit  length 


p^Pt 

m 

m 


coefficients  associated  with  piston  theory  aerodynamics 
3x3  constitutive  matrix  for  a  layer 
material  properties  matrix 


Qi  generalized  force 

{Q}  vector  of  generalized  forces 

{Qdamp]  aerodynamic  stifftiess  and  damping  matrices 
qp,{(l}  generalized  displacement  and  the  vector  of  generalized  displacements  for  panel 


qjlutU  qcrit 

R,S 


r. 


u 


V 

W'(x,y) 

w 

qw,pw 

[a] 


5w 


0 


^flutu  ^crit 

e 

P<^  y  pm 


CO 


¥ 


flutter  dynamic  pressures  with  and  without  damping,  respectively 
coefficients  of  front  line  or  aft  line  of  panel 
thickness  of  layer  i 

kinetic  energy  of  the  system 

coefficient]  in  the  polynomial  series  for  layer  i 

total  potential  energy  of  the  system 
shape  dependent  coefficients  of  Fb 

free  stream  velocity 

fnst  order  state  space  system  matrices 

potential  energy  of  the  panel  due  to  in-plane  loads 
panel  elastic  out  of  plane  displacement 
coefficients  of  polynomial  elements  of  [a] 

matrix  of  function  derivatives 

virtual  work 

yaw  angle  of  the  flow 

normalized  flutter  dynamic  pressures  with  and  without  damping 
eigenvalue  of  the  aeroelastic  system 
fiber  orientation  angle 

free  stream  density  and  panel  material  density 
frequency 

right  and  left  eigenvectors  of  the  generalized  eigenvalue  problem 
panel  skew  angle 


V 


normalized  flutter  frequency 


^crit 

subscripts 
A  -  aft  (rear) 
F  -  front. 

L  -  left 
R  -  right. 


VI 


ACKNOWLEDGMENTS 


Firstly  I  thaiik  my  Lord  Jesus,  who  has  blessed  me  with  die  ability  and  opportunity  to  complete 
this  degree.  I  would  like  to  express  my  thanks  to  Professor  Eli  Livne  for  his  tremendous  vision,  guidance, 
and  support  during  the  length  of  my  graduate  studies  and  the  preparation  of  this  work.  His  vast 
knowledge  in  all  areas  of  engineering  has  left  a  lasting  mark  on  my  education  and  is  greatly  appreciated. 
To  my  wife  Amy-thank  you  so  much  for  your  support  and  love  over  the  last  several  months.  You  kept  my 
sanity  intact. 


CHAPTER  1 


INTRODUCTION 


A  substantial  amount  of  experience  and  knowledge  has  been  accumulated  over  the  last  forty  years 
regarding  the  aeroelasticity  of  panels  in  supersonic  flow  (Refs.  1-4).  Different  numerical  solution 
techniques  have  been  used,  including  exact,  as  well  as  approximate  Galerkin,  Rayleigh-Ritz,  and  Finite 
Element  techniques  (Refs.  5-7).  The  problem  has  practical  implications  associated  with  the  design  of  high 
speed  aerospace  vehicles.  As  an  aeroelastic  research  problem,  it  is  rich  and  interesting,  encompassing  a 
wealth  of  phenomena.  These  include  linear  and  nonlinear  stability  and  dynamic  response,  dynamics  of 
systems  with  random  parameters  (Ref.  8),  and  interactions  between  static  and  dynamic  instabilities  in  the 
presence  of  in-plane  loads  and  thermal  effects  (Refs.  9-10).  Panel  flutter  has  been  used  to  study 
applications  of  composite  materials  (Refs.  11-12),  transverse  shear  effects  (Refs.  13-14),  active  aeroelastic 
control  using  strain  actuators  (Refs.  15-16),  chaotic  dynamics  (Refs.  17-18) ,  and  order  reduction  in 
unsteady  aerodynamics  (Ref.  19),  Most  of  the  studies  in  the  vast  panel  flutter  literature,  however,  are 
confined  to  rectangular  panels.  Solution  techniques  and  results  for  quadrilateral  and  trapezoidal  panels 
have  been  discussed  in  only  a  small  number  of  articles  (Refs.  20-23),  Skin  panels  in  typical  aerospace 
structures,  however,  are  very  often  trapezoidal  in  shape.  Moreover,  in  the  course  of  optimization  of  such 
aerospace  structures,  internal  ribs,  spars  and  stiffeners  may  be  moved  to  form  trapezoidal  skin  panels,  and 
these  panels  may  change  shape  during  optimization  in  addition  to  changing  material  properties  and 
thicknesses.  The  capability  to  efficiently  evaluate  the  aeroelastic  stability  of  trapezoidal  skin  panels  under 
combined  in-plane  loads,  as  well  as  sensitivities  with  respect  to  sizing,  material  and  shape  design 
variables,  constitutes  an  important  building  block  in  any  overall  structural/aeroelastic  optimization 
capability.  For  optimization  strategies  based  on  non-linear  programming  and  approximation  concepts 
(Refs.  24),  evaluation  of  alternative  approximation  techniques  for  aeroelastic  constraints  is  required,  since 
very  little  experience  exists  in  this  area. 

The  optimization  of  panels,  subject  to  aeroelastic  constraints,  has  been  studied  in  cases  involving 
isotropic  and  composite  construction  (Refs.  25-29).  These  studies  are  usually  limited  to  the  treatment  of  an 
isolated  panel,  excluding  its  interaction  with  the  structure  containing  it.  In  our  effort  to  develop  effective 
aeroservoelastic  synthesis  techniques  for  stressed  skin  aerospace  structures,  we  have  ventured  into  the  area 
of  airframe  shape  optimization  in  order  to  make  rigorous  design  optimization  available  to  the  designer  at 
an  early  stage  of  the  design  process,  where  overall  shape  of  the  vehicle  is  still  evolving.  Analytic  shape 
sensitivities  and  approximations  have  been  developed  for  wing  box  structural  modeling,  integrated  wing 
box  /  panel  buckling  analysis,  and  unsteady  wing  aerodynamics  in  both  subsonic  and  supersonic  flight 
(Refs.  35-38).  The  present  work  will  focus  on  the  integrated  wing  box  /  panel  flutter  analysis  and 


sensitivity  problem  with  an  emphasis  on  the  needs  of  planform  shape  optimization.  This  includes  the 
sensitivity  analysis  of  panel  flutter  in  the  case  of  the  isolated  panel,  subject  to  given  in-plane  loads. 

The  thesis  opens  in  Chapter  2  with  a  derivation  of  the  equations  of  motions  for  the  panel  flutter 
problem  based  on  a  Ritz  approximation.  Chapter  3  discusses  the  panel  modeling  and  the  simple 
polynomials  used  as  admissible  functions  in  the  Ritz  analysis.  Chapters  4  through  7  develop  the  mass, 
stiffness,  geometric  stifftiess,  and  aerodynamic  matrices.  Analytic  sensitivities  with  respect  to  panel 
shape,  thickness,  and  fiber  orientation  are  derived  in  Chapter  8.  Chapter  9  presents  the  verification  and 
results  of  the  present  technique  and  Chapter  10  provides  the  conclusions  of  the  work. 


CHAPTER  2 


PANEL  FLUTTER  EQUATIONS  OF  MOTION 

2.1  Introduction 


The  panel  flutter  equations  and  stability  analysis  are  developed  in  this  chapter.  The  equations  of 
motion  are  derived  using  energy  methods  and  Ritz  approximation.  The  formulation  uses  the  assumptions 
of  classical  plate  theory  for  symmetric  composite  laminates.  A  State  Space  modeling  approach  is  used  to 
develop  the  stability  analysis  of  the  system,  and  the  flutter  problem  is  shown  to  be  a  linear  generalized 
eigenvalue  problem.  The  determination  of  critical  dynamic  pressures  calculated  with  and  without 
aerodynamic  damping  is  discussed  as  well  as  the  different  types  of  aeroelastic  instability. 


2.2  The  Lagrange  Equations 


The  standard  Lagrange  equations,  for  a  system  with  N  degrees  of  freedom,  is  given  as: 


d  t 


( 


d  T 


d  (T-U) 
a  qj 


=  Qj  j  =  l..N 


(2.1) 


where  qj  are  the  generalized  coordinates,  U  is  the  total  potential  energy  of  the  system,  T  is  the  kinetic 
energy,  and  Qj  (the  generalized  forces)  are  defined  by  the  virtual  work  of  non-conservative  forces: 


5  w  =  ^Qj  5  qj 


(2.2) 


In  plate  theory,  the  potential  energy  is  usually  derived  as  the  sum  of  the  strain  energy  due  to 
bending,  U,  and  the  potential  energy  due  to  in-plane  loads,  V.  Thus  Equation  2.1  can  be  written  as 


a  (u+v-T) 

d  qj  d  t 


^  d  T 


=  Qj  7  =  1...A^ 


(2.3) 


For  a  symmetric  laminate  the  strain  energy  due  to  bending  (Ref.  52)  can  be  written  in  the  form 


4 


1^4  21^4  }d 

where  (x^y^t)  is  the  vertical  displacement  of  the  panel  and  [D]  is  the  symmetric  3X3  bending  stiffness 
matrix  defined  by 

D^=\z'Qydz  (2.5) 

z 


yy 


2W‘ 


dxdy 


(2.4) 


The  potential  energy  of  in-plane  loads,  Nj,  due  to  the  deflection  W‘(;i,y,t)  is  (Ref.  52) 


Finally,  the  kinetic  energy  of  the  panel  is  given  by 


where  pm  is  the  material  density  per  unit  volume  and  h(x,y)  is  the  panel  thickness. 

2.3  Ritz  Formulation 


(2.6) 


(2.7) 


The  unknown  elastic  deformation  W\x,y,t)  is  approximated  by  a  series  of  admissible  functions; 

N 

W\x,y,t)  =  Yjfi(x,y)^M')  (2.8) 

The  functions  f,  (x,y)  satisfy  the  geometric  boundary  conditions  and  qi(t)  are  the  generalized  displacements. 
In  matrix  form  the  deformation  is: 
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'*'i,,,0={A<x,y)  /2(x,y) 


^l(0 


Q  yV  (0 


(2.9) 


The  first  derivatives  can  be  expressed  as: 


K-l. 

'f\.x 

fl,  . 

> 
_ 1 

K.J- 

h.,  ■ 

- 1 

^1 

^2 

J 


(2,10) 


and  the  second  derivatives  are  given  by: 


W\xx  ' 

^  _ 

fljcx  flja 

fl.yy  fl.yy 

•  /v.xx 

*  f N,yy 

^1 

^2 

■  =  [''3X9} 

2W\xy 

}fuxy  2/,,^  . 

••  ^fN,xy^ 

(2.11) 


These  Ritz  approximations,  used  to  reduce  the  order  of  the  system  from  infinite  to  N,  can  now  be 
substituted  into  the  Lagrange  equations. 


2.4  Strain  Energy 


The  partial  derivative  of  the  strain  energy  (Eq.  2.4)  with  respect  to  the  generalized  coordinate  q. 


is: 
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W', 

1W,\ 


1 

yy 


xy 


+[W'.L 


(Twi 

rr  , 


a 


w. 


2W  ‘ 


(2.12) 


\ixdy 
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From  equation  2.8  the  second  derivatives  of  can  be  given  as: 


H'.L  =  I  =  I  W.;,  =  S  /,,9i 


(2.13) 


from  which  follows 


ajC 

3^, 


d  qj 


d  W, 


a  qj 


—  f 

J  J.xy 


(2.14) 


substituting  Eq.  2.14  into  Eq.  2.12  gives: 


a  U 
a  qj 


4JJ[4 


yy 


'  w  ‘ 

’XX 

J 

W  ‘ 

J  ’yy 

+  K  n 

_2iy4_ 

(2.15) 


Notice  that  the  quantity 


yy 


2W4Jd] 


fj.X. 

fj.yy 

y,. 


is  a  1x1  matrix  and  is  therefore 


symmetric.  Thus 


\  fj,xx 

fj,xx 

[w4  ix,;,  2w4);d1 

[IV.L  IV,;,  2W,'^Id] 

fj.yy 

W'.L 

W' 


Substituting  Eq.  2.16  simplifies  2.15  to: 


(2.16) 


a  u 


h. 


2/,H 


yy 


2W,1 


\dxdy 


(2.17) 


For  all  j=l...N  Eq.  2.17  becomes: 
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fl,yy 
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2/v.jo’. 
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dxdy 


a  Qn 

Substituting  Eq.  2.11  gives: 

fa  f/l 


a  qi 

a  u 
a  ^2 

a  {/ 


a  q^ 


2.5  Potential  Energy  Due  to  In-plane  Forces 
The  partial  derivative  of  V  (Eq.  2.6)  with  respect  to  the  generalized  coordinate  qj  is: 


a  qj 

where  M= 


dxdy 


(2.18) 


(2.19) 


(2.20) 


From  equation  2.8  the  first  derivatives  of  (x,y,t)  can  be  given  as: 


(2.21) 
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from  which  follows 

d  W  ‘ 


=  fj. 


d  w,; 

d  qj 


—  f  ■ 

J  j.y 


(2.22) 


substituting  Eq.  2.22  into  2.20  gives: 


(2.23) 


As  in  the  derivation  of  the  strain  energy,  the  second  term  in  the  integral  is  symmetric  because  it  is  a  1x1 
matrix.  Therefore  Eq.  2.23  can  be  written  as: 


(2.24) 


Forallj=l...N 


(2.25) 


Finally  substituting  Eq.  2.10  into  Eq.  2.25  gives: 
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idV] 
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d  V 


P  ^2  ,1 


d  V 
5  Qn  _ 


(2.26) 


2.6  Kinetic  Energy 


The  kinetic  energy  in  this  problem  (Eq.  2.7)  is  a  function  of  the  time  derivatives  of  the 

a  T 

generalized  coordinates  only,  therefore - =  0 . 

a 

The  derivative  of  W\x,y,t)  with  respect  to  time  can  be  written  as: 


a  w' 

d  t 


(2.27) 


Substituting  this  expression  into  Eq.  2.7  results  in: 


dxdy 


-<=i  J 


Differentiating  with  respect  to  gives: 


a  T 
^  <!],• 


N  \ 


dxdy 


.  1=1  y 


(2.28) 


(2.29) 


Differentiating  with  respect  to  time  gives: 


21  Virtual  Work  Due  to  Transverse  Loads 

The  virtual  work,  5w  done  by  a  transverse  force,  ,  on  the  panel  is  the  force  over  the  panel 
area  multiplied  by  the  virtual  distance: 

5  w  =  jj5  q^^y^dxdy  (2.32) 

Using  Eq.  2,8  the  virtual  displacement  can  be  expressed  as: 

5H'‘=|;/,5  4, 

i=l 


Combining  the  above  equations  gives: 


(2.33) 


5  vt;  =  jJX/,5 

/=1 

By  the  definition  of  Qi  (Eq.  2.2),  the  vector  {Q}  can  be  expressed  as: 


te}=|J 


7i 

fN 


‘l(.,y)dxdy  =  q^,^y)dxdy 


(2.34) 


^(x  y)  transverse  load,  but  in  the  context  of  panel  flutter  it  is  the  pressure  applied  to 

the  plate  by  the  aerodynamic  flow.  Therefore  the  aerodynamic  contribution  to  the  system  will  be  solely 
through  the  {Q}  vector. 


2.8  Equations  of  Motion 


Recall  that  the  Lagrange  equation  was  given  by: 


d 


{U+V) 

dq^ 


d 

+ - 

8  t 


^  d  t'^ 
7  J 


=  Qi 


Substituting  the  derived  equations  above  gives; 


jj  [^3  f  [^p3  ]{q}dxdy  +  JJ  [F,  f  [)V  JF^  ]{^  }dxdy  +  JJ  p  „;i[Fi  f  [Fj 


Because  the  generalized  coordinates,  {q},  are  not  functions  of  x  and  y  they  can  be  placed  outside  of  the 
integrals.  The  equations  of  motion  can  then  be  expressed  as: 


(2.36) 
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where  [M]  is  the  NxN  mass  matrix: 


[K]  is  ttie  NxN  bending  stiffness  matrix: 

[K\^\\[F,nDlF,Wy 


(2.37) 


(2.38) 


and  [Ko]  is  the  NxN  geometric  stiffness  matrix: 

[Koh\\[F^^[N\F^}lxdy  (2.39) 

2.9  Generalized  Eigenvalue  Problem 

It  win  be  shown  in  Chapter  7  that  when  piston  theory  aerodynamics  are  used  for  the  transverse 
load  contribution,  the  nonconservative  aerodynamic  forces  are  expressed  as: 

{Q}= K^} 


Substituting  this  into  Eq.  2.36  gives  the  equations  of  motion  as: 

]|?|  =  joj 

or 


(2.41) 


(2.42) 


It  should  be  noted  that  the  model  presented  above  does  not  include  any  structural  damping.  This  can  be 
added  in  the  form  of  viscous  damping  in  the  |cj matrix  (Ref.  31). 
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At  a  given  altitude,  Mach  number,  and  corresponding  speed,  and  for  a  given  set  of  in-plane  loads 
obtained  from  the  wing  box  solution  (Ref.  37),  the  poles  of  the  system  can  be  found  to  determine  if  the 
panel  is  stable  or  unstable.  Applying  the  Laplace  transform  to  Eq.  2.42  gives; 

5"  M  +  5C  -f  }  =  {0}  (2.43) 


This  is  a  quadratic  eigenvalue  problem  which  can  be  solved  by  converting  it  to  a  first  order  problem. 
Let  -[x;  (5)}  =  {^(5)}  and  {^2  (i')}  =  i{^(5')}then  Eq.  2.43  becomes: 


5{x,}={x2} 

^  |^{^2  }  =  -  ^]{^i }  -  p]{^2  } 


(2.44) 


or 


(2.45) 


If  we  have  N  generalized  coordinates,  {q},  then  M  ,  cj,  and  j^j 


are  all  NxN  matrices  and  the  first 


order  equations  are  a  2Nx2N  generalized  linear  eigenvalue  problem  of  the  form: 


(2.46) 


where  X  are  the  eigenvalues  (poles)  of  the  system  and  {([)}  are  the  corresponding  eigenvectors.  The 
eigenvalues  and  eigenvectors  are  found  using  the  QZ  algorithm  (Ref.  42). 

2.10  Stability  Boundaries 


At  the  given  flight  conditions  and  in-plane  loads,  a  panel  is  stable  if  all  of  its  poles  (eigenvalues) 
reside  in  the  left  side  of  the  Laplace  plane.  Each  pole  contributes  to  a  transient  motion  in  the  panel  of  the 
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form  .  If  a  pole  is  written  as  ?i=a+jco,  then  it  contributes  a  motion  of  the  form  =  e^^^Ccoscot  + 
jsincot).  A  pole’s  contribution  to  the  panel  displacement  grows  with  time  if  its  real  part,  a,  is  greater  than 
zero  and  therefore  causes  instability  in  the  panel.  A  complex  pole  with  a  positive  real  part  represents  an 
oscillatory  instability,  or  flutter.  A  purely  real  pole  that  is  greater  than  zero  is  non-oscillatory  and 
corresponds  to  divergence  or  buckling.  Thus,  eigenvalue  analysis  can  capture  instability  due  to  either 
buckling  or  flutter. 

The  point  at  which  a  panel  becomes  unstable  is  usually  found  by  assuming  a  given  Mach  number 
and  set  of  in-plane  loads.  The  dynamic  pressure  is,  then,  slowly  increased,  keeping  all  other  parameters 
constant.  Eigenvalue  analysis  is  repeated  for  increasing  dynamic  pressures  until  a  pole  crosses  into  the 
right  hand  side  of  the  s-plane  (Fig.  2.1).  This  is  equivalent  to  flying  at  a  constant  Mach  number  while 
decreasing  altitude.  The  dynamic  pressure  at  which  the  panel  flutters  is  labeled  qnutter- 


In  the  panel  flutter  literature,  it  has  been  observed  that  due  to  limitations  of  Piston  Theory,  the 
associated  aerodynamic  damping  can  sometimes  lead  to  counter  intuitive  /  unreliable  results.  Structural 
damping  effects  are  also  quite  hard  to  model  accurately ( Ref.  31).  For  a  conservative  stability  boundary, 

then,  the  structural  and  aerodynamic  damping  are  often  ignored  by  setting  |c j  =  [O]  in  Eq.  2.42.  This 

results  in  a  model  with  no  damping,  and  the  stability  analysis  has  to  be  modified  accordingly. 


As  long  as  the  system  is  stable  its  poles  are  purely  imaginary.  Only  after  flutter  or  buckling 
occurs  does  any  root  have  a  real  part  (Fig.  2.2).  In  the  case  of  flutter,  two  poles  that  start  out  as  natural 
frequencies  will  coalesce  as  the  dynamic  pressure  is  increased  to  the  point  where  they  become  complex 
conjugates  of  the  form  ±a  +  jco  (Fig.  2.2a).  The  corresponding  dynamic  pressure  is  caUed  qcriticai- 
Buckling  or  divergence  occurs  when  one  root  becomes  a  positive,  non-oscillatory  root  (Fig.  2.2b).  The 
present  analysis  (and  associated  computer  program)  predict  both  qcriticai  and  qnutter- 


Figure  2.2  Stability  Analysis  in  the  Absence  of  Damping 


CHAPTERS 


POLYNOMIAL  MODELING 

3.1  Introduction 

Simple  polynomials  for  modeling  and  Ritz  approximation  basis  functions  are  discussed  and 
developed  in  this  chapter.  The  motivation  for  using  simple  polynomials  is  discussed  followed  by  a 
presentation  of  the  wing  box  and  panel  models  used  in  the  present  analysis.  Finally,  admissible 
polynomials  for  simply  supported  panels  are  developed. 

3.2  Motivation  for  Simple  Polynomials 

Admissible  functions  based  on  simple  polynomials,  as  is  well  known,  lead  to  ill-conditioned 
system  matrices  (Refs.  32-34).  As  higher  order  polynomials  are  introduced  to  better  approximate  the 
solution,  high  and  low  order  terms  appear  simultaneously  in  the  system  matrices,  leading  to  a  large 
difference  in  orders  of  magnitude  of  matrix  terms.  On  finite  word-length  computers  this  leads  to  ill- 
conditioned  linear  equation  and  eigenvalue  solutions.  If  convergence  to  the  solution  is  not  obtained  with 
low  order  polynomials,  the  addition  of  higher  orders  will  lead  to  ill-conditioning  and  the  solution 
technique  will  likely  fail. 

However,  there  are  several  advantages  to  simple  polynomials  if  ill-conditioning  can  be  avoided. 
Simple  polynomials  allow  for  integrals  in  the  system  matrices  to  be  performed  analytically,  making 
numerical  quadrature  unnecessary.  This  greatly  reduces  the  computation  time  required  for  solution. 
Simple  polynomials  have  been  successfully  used  as  admissible  functions  in  plate  vibration  and  shells.  In 
the  context  of  wing  optimization,  simple  polynomials  have  been  used  for  deformation,  stress,  and  mode 
shapes  in  wings  of  complex  planforms  (Ref.  35),  Finally,  Ref.  37  shows  that  panel  buckling  analysis  can 
be  successfully  integrated  with  wing  box  structural  analysis  using  simple  polynomials. 

The  greatest  advantage  of  simple  polynomials  for  the  present  problem  of  panel  flutter  is  the 
ability  to  perform  analytical  sensitivities.  Refs.  35  and  37  show  that  it  is  possible  to  easily  obtain  very 
accurate  closed  form  sensitivities.  In  design  oriented  structural  analysis  there  is  a  need  to  analyze  the 
flutter  stability  boundary  of  the  wing  panels  when  their  shapes  and  other  design  variables  are  changing. 
Analytical  sensitivities  for  the  stability  of  wing  panels  that  can  be  provided  quickly  to  a  designer  provides 
an  important  link  in  wing  and  airframe  optimization. 

The  success  of  previous  experiences  using  simple  polynomials  for  analysis  and  analytical 
sensitivities  combined  with  the  computational  speed  such  formulation  allows  is  the  basis  for  using  simple 
polynomials  in  the  present  panel  flutter  analysis. 
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3.3  Modeling 

In  the  polynomial  based  wing  box  analysis,  the  thickness  of  layers  of  fibers  in  different  directions 
is  given  by  simple  polynomials.  For  layer  i  out  of  Nt  layers  in  a  panel 

,i=l,NL  (3.1) 

k  ^  1 

The  powers  and  are  x  and  y  powers  of  the  term  of  the  thickness  series  for  the  i* 

layer.  The  coefficients  T^  serve  as  sizing  type  design  variables.  Unlike  many  studies,  in  which  wing 

trapezoidal  segments  are  transformed  into  a  unit  square  for  numerical  analysis,  here  the  thickness 
polynomial  is  given  in  terms  of  the  physical  x,y  coordinates  (Fig.  3.1). 


X 


Air  Flow 


Figure  3 . 1  W ing  Planf orm 


The  overall  thickness  of  skin  panels  is  the  summation  of  the  Nl  layers: 

V, 


L  i  (mh  {nh 


(3.2) 


In  a  similar  manner,  cap  areas  for  the  spars  and  ribs  of  the  wing  box  model  are  also  expressed  as 
simple  polynomials  of  either  x  (for  ribs)  or  y  (for  spars)  (Refs.  32-34).  Depth  of  the  wing  box  is  also  given 
by  a  simple  polynomial  in  x  and  y,  to  be  defined  later. 


3  A  Admissible  Functions 
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An  admissible  function  will  be  any  polynomial  we  choose  that  satisfies  the  simply  supported 
boundary  conditions  on  the  perimeter  of  the  panel.  Figure  3.2  shows  a  trapezoidal  panel  defined  by 
coordinates  of  its  vertices  in  the  x,y  axes.  The  subscripts  L  and  R  denote  left  and  right  sides,  respectively. 
Tlie  subscripts  F  and  A  denote  front  and  aft  lines,  respectively,  and  xp  and  xa  are  the  x  coordinates  of  the 
forward  and  rear  points  on  a  line  parallel  to  the  sides  of  the  panel. 


Based  on  Fig.  3.2,  we  can  write  the  following  equation  for  points  on  the  front  line: 

^  (y)  =  ^IfiZrJLeeIl,^ + (If/?  =  R  +s  y 

^  V  — V  V—  a;  r  r 


^R-yi 


yR-yi 


(3.3) 


In  a  similar  way,  on  the  rear  line,  expressing  xa  in  terms  of  y  leads  to: 


yR 


(3.4) 
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The  function  h  {x, y)  =  ^x-xp ()^) Jx - x^  ( j)|y  -  Jt  -  y/j  ]  satisfies  the  zero 

displacement  boundary  conditions  on  the  circumference  of  the  panel.  Using  Eqs.  3.3  and  3.4,  and 
expanding  in  terms  of  x  and  y  yields 

'^S(x,y)  =(f'l  +V^’c  +  V^y*V^x^-.V^xy-^V^yh 


or 


3  nf<  6  tn^:  n)  3  6  nf^+n) 

(x.y)  '  I,'';-  ^  =, 111, ^ 


(3.6) 


where  the  constants  U  and  V  are  given  in  terms  of  panel  vertex  coordinates  by: 


K,  =  R^Rp.V^  =  -(r^  +  Rf  ).V3  =  {r^S^^R^S^) 

=  l,v^^  = -(5^  +  5^),y^  =  5^5^ 

Powers  of  x  and  y  corresponding  to  constants  U  and  Vare  given  in  Tables  3.1  and  3.2. 


Table  3.1 _ Constants  Ui  and  their  corresponding  power  of  y. 


i 

Ui 

u 

“i 

1 

Ui 

0 

2 

U2 

1 

3 

U3 

2 

Table  3.2 _ Constants  Vj  and  their  corresponding  powers  of  x  and  y. 


j 

Vj 

V 

mj 

V 

nj 

1 

0 

0 

2 

V2 

1 

0 

3 

V3 

0 

1 

4 

V4 

2 

0 

5 

V5 

1 

1 

6 

V6 

0 

2 
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Any  polynomial  multiplying  Fb  will  be  an  admissible  function  satisfying  the  boundary  conditions 
of  a  simply  supported  panel.  Therefore  the  transverse  displacement  of  the  panel  can  be  written  as: 

1  %  mW 

where  the  coefficients  q  p  are  the  generalized  displacements.  Substituting  the  Eq.  3.6  for  Fb  we  can 


write: 


rl, 


N, 


W 


3  6 


p=\  i=iy=i 


V,V.x 


(mV+m^)  (nW+ny+n^) 


(3.9) 


The  admissible  functions  are  thus  expressed  in  terms  of  simple  polynomials,  where  the  p^ 
admissible  function  (Eq.  2.8)  is  given  by: 


(x,y) 


3  6  (mV+m^)  (nf+n^+n^) 

=  I  I  UjV-x  J  ^  y  '  J  ^ 
i=lj=l  '  J 


(3.10) 


CHAPTER  4 


MASS  MATRIX 

4.1  Introduction 

Equations  for  the  terms  of  the  mass  matrix  are  derived  in  this  chapter.  It  is  shov^n  that  the  terms  can 
be  expressed  as  linear  combinations  of  area  integrals  of  polynomials. 

4.2  Mass  Matrix  Terms 
The  mass  matrix  was  derived  in  Chapter  2  as: 

(4,1) 

area 

where  pm  is  the  material  density  per  unit  volume,  h(x,y)  is  the  panel  thickness,  and  the  matrix  [Fi]  contains 
the  admissible  functions: 

['",]=[/,  A  /»1  (4.2) 

Substituting  the  admissible  functions  into  [M]  gives  the  individual  terms  as: 

=  \jpmhfrfsdxdy  (4.3) 

area 

This  equation  shows  that  the  mass  matrix  is  symmetric. 

The  panel  thickness,  h,  is  expressed  as  a  summation  of  the  individual  layers’  thickness.  As  shovm  in 
Chapter  3  the  thickness  is: 

1=1  i=i 

where  Nl  is  the  number  of  layers  in  the  panel  and  Ni  is  the  number  of  thickness  terms  in  each  layer.  If  the 
thickness  of  all  layers  are  expressed  by  the  same  order  of  complete  polynomial  then  all  layers  will  have  the 

same  number  of  terms,  Nt ,  and  the  same  powers  •  This  restriction  simplifies  the  formulation 

of  the  panel  thickness  while  not  taking  away  from  the  flexibility  for  representing  complex  panel  thickness. 
A  new  vector,  7^  is  defined  to  add  all  of  the  terms  associated  with  the  same  powers  of  x  and  y: 

%=^n  (4,5) 

1=1 

Now  the  total  panel  thickness  is  simplified  as: 


nx” 


(4.4) 


L 


22 


jt=i 


Substituting  this  into  equation  (4.3)  leads  to 

^rs=9miT^\\  dxdy 


(4.6) 


(4.7) 


area 


Recall  that  the  p  admissible  function  was  given  in  Chapter  3  in  terms  of  the  coefficients  Ui  and  Vj  as 


3  6 


a:'"-’*"'’ 


(4.8) 


/=1  ;=l 


Substituting  into  equation  (4.7)  leads  to  the  explicit  form  of  the  mass  matrix  terms: 

=  p.liiii  Wjt'.n  }f^"y'‘ixdy 


(4.9) 


k=l  i=l  J=l  «=1  ;;=1 


where 


m  =  m''j  +  nijj  +  +  m[ 

n  =  n“  +n"  +n]  +n]j  +<  +<  +n[ 

Note  that  all  of  the  elements  of  the  mass  matrix  are  a  linear  combination  of  area  integrals  of  simple 
polynomials 

=  \\x"'y"(bcdy  (4.10) 

area 

which,  for  trapezoidal  panels,  can  be  carried  out  analytically,  as  shown  in  detail  in  Ref  .35.  It  will  be 
shown  that  all  of  the  system  matrices  will  be  combinations  of  the  same  family  of  integrals  Itrcidji)  • 

The  mass  matrix  is  dependent  only  on  the  thickness  terms  through  7]^  and  the  panel  shape 

variables  through  Uj  and  Vj.  The  area  integrals  depend  only  on  the  panel  shape  variables  through  the 
limits  of  integration  as  shown  in  Ref  35. 


CHAPTER  5 


STIFFNESS  MATRIX 


5  A  Introduction 

Expressions  for  the  terms  of  the  stiffness  matrix  [K]  are  derived  in  this  chapter.  It  is  shown  that 
the  stiffness  terms  can  be  expressed  as  combinations  of  integrals  of  simple  polynomials.  Quasi 
homogeneous  or  orthotropic  construction  is  assumed. 


5.2  The  In-plane  Stiffness  Matrix 

As  shown  in  Chapter  2,  the  equation  for  the  Ritz  formulation  of  the  stiffness  matrix  is: 

<51> 

area 

where  [F3]  contains  the  second  derivatives  of  the  admissible  functions  and  [D]  is  the  bending  stiffness 
matrix.  Expressions  for  the  in-plane  stiffness  matrix,  [A],  the  bending  stiffness  matrix  and  the  matrix 
[F3]  are  required  for  an  expression  for  [K]. 

The  in-plane  stiffness  matrix  is  defined  as 

[^]  =  J  [C(z)yz  (5.2) 

Z 

where  [Q(z)]  is  the  material  property  matrix.  [Q(z)]  can  be  expressed  in  terms  of  material  invarients  and 
the  layer  fiber  orientation  angle  by  defining  [Qo]  through  [Q4]  (Refs.  37,51)  as: 


[£’]= 


'U, 

U, 

o' 

'1 

0 

o' 

'  1 

-1 

o' 

0 

+u. 

0 

-1 

0 

cos(20)  +  1/3 

-1 

1 

0 

_  0 

0 

^5. 

_0 

0 

0 

_0 

0 

-1 

u. 


'0 

0 

r 

'0 

0 

n 

0 

0 

1 

sin(20)  +  t/3 

0 

0  -1 

_1 

1 

0_ 

_1 

-1 

oj 

sin(46) 


cos(40)  + 


(5.3) 


The  polynomial  description  of  panel  thicknesses  in  terms  of  the  global  x-y  coordinates  was  given  in 
Chapter  3  as  a  summation  of  the  thickness  terms  of  Nl  layers: 

/=i  i=i  k=i 


Nr  N, 


YYnx’ 


=  t 


l(x,y) 


Ni{x,y) 


L 

='Lh 


(x.y) 


(5.4) 
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For  a  panel  containing  Nl  layers  of  fibers,  the  in-plane  stiffness  matrix  [A]  can  be  expressed  in  terms  of 
individual  layer  thicknesses  and  fiber  orientation  angles  by  substituting  Eqs.  5.3  and  5.4  into  Eq.  5.2: 


Nl 

[A]=  .lJ[(2o]+[!2i]cos2ej-  +[e2]cos4e^-  +[23]sin2e-  -b[!24]sin4ej]  f. 


(5.5) 


Let  a  material  and  fiber  orientation  dependent  matrix,  [eo  j)j,  a 3x3 

[5(6/)]  =  [[2o]+[2i]cos2e-  -i-[22]cos4e-  -i-[e3]sin2ej-  +[24]sin4ej] 


3x3  matrix,  be  defined  as: 


(5.6) 


The  in-plane  stiffness  matrix  [A]  can  now  be  expressed  in  terms  of  the  sizing  design  variables,  T^j ,  and 
fiber  directions  as  a  polynomial: 


[x]=  s  i  [ete,)’  ■ 

i=U=lL  '  J 


(5.7) 


5.3  The  Bending  Stiffiiess  Matrix 

For  unidirectional,  orthotropic  or  quasi  homogeneous  laminates  the  in-plane  and  bending 
stiffness  matrices  are  related  through  (Ref.  54) 

.2 


[D]  =  [A] 


hL 

12 


(5.8) 


Using  Eq  5.4  to  express  h^  in  terms  of  thickness  design  variables,  double  summation  is  needed. 
The  indices  li  and  I2  are  used  for  summation  of  polynomial  terms  associated  with  each  layer,  as  follows; 


N, 


\x,y) 


=  .S.  I  r; 


/1=1/1=1 


l\ 


(m^l) 


iV, 


^•2 


=  X.  X  T. 


il 


/2=1/2=1 


12 


(5.9) 


The  bending  stiffness  matrix,  [D],  can  now  be  written  in  polynomial  form  as 
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[d]=[4 


12 


—  I  0(0,)  I  I  I  I  I  r' 
12t=lL  '  -tl=lj2=U=in=l/2=l  ^ 


-pi\  rpi2 

'Hi 


ini'  ^rnli^+niji^  ( n{i  +/2f5^ 
k  l\  12  \  k  l\  12 

y- 


(5.10) 


The  dependence  of  [D]  on  the  sizing  design  variables  (thickness  coefficients)  and  fiber  angles  is 
now  expressed  in  explicit  form. 


5.4  The  fFjl  Matrix 

Recall  from  Chapter  2  that  the  [F3]  matrix  is  defined  as: 


fl.xx 

fl^xx 

In^xx 

[F,]  = 

f l,yy 

f 2,yy 

f N  ,yy 

(5.11) 

2/i.^ 

2/2,^  * 

••  '^fn^xy 

The  p  admissible  function  is  given  as: 

3  6 


pi^^y) 


1=1  ;=i 


(5.12) 


Second  derivatives  of  the  fp  functions  are: 


3  6  V  V  w  K+<-2)  (/f+n^+n^) 

fp,xx=X,J.ViVji.m'’:+nX)(.m)+nX-\)x  J  ^  y‘  ^  ^ 

1=1 7=1  J  J  J 


(5.13) 


Note  that  fp^x  =  0,  if  m  ^  +  m  ^  <  1  •  (representing  second  derivatives  of  zero  order  or  first  order 
terms)  Similarly 


(nf +n;;+«^)(nf +n}+n 


^  {my+mp  (npnyn^-2) 
p  y 


(5.14) 
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where  fp^yy  =0,  if  +n^j +np  <l .  Finally 


'’+m^)(nf+n'’j+n^)x 


(m^+m^-1)  (n^+n^-+n^-i) 


(5.15) 


The  term  fp,xy  is  set  to  zero  in  two  cases: 

2  fp^xy  -  0’  if  or  +n^  +  tip  =  0 


The  q,p  element  of  [F3]  can  be  written  as: 
f  _  y  y  F^r  ”  v 


(5.16) 


where  the  coefficients  and  corresponding  powers  of  x  and  y,  m/,f^  and  nf:^  respectively,  are 

V  V  y 


shown  in  Table  5.1. 


Table  5.1  Coefficients  and  powers  of  polynomial  tenns  in  the  [F3]  Matrix 

(If  any  of  the  power  of  x  and  y  from  columns  three  and  four  is  less  then  zero,  the  corresponding  F  ^ 

element  is  set  to  zero) 


q 

(Row  of  [F3]  Matrix) 

pqp 

ij 

mf9P 

nfSP 

(Coefficient) 

(Power  of  x) 

(Power  of  y) 

1 

U  iVj(m  J  +mp)- 

(mj^+m^-l) 

mj  +m^-2 

a■^  +nj  +np 

2 

UiVj(n|*  +  nJ+np- 

(nH  +  ny  +  nJ^-1) 

mj  +mp 

n^  H-  Oj  +np  -2 

3 

u,y(ra;  +  m;). 

mJ 

_Lr.  V  1  „w  . 

n|  +nj  H-Dp  -1 

(n“  +ny  +n*) 

V  '  J  P  / 

5.5  The  Stiffness  Matrix 
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Equation  5.16  together  with  Eq.  5.10  are  substituted  into  the  equation  for  the  stiffness  matrix 
(Eq.  5.2).  The  r,s  term  of  the  stiffness  matrix  is: 

"ai  bsH>’  <3  ”» 

and  the  final  expression  for  the  Krs  element  in  polynomial  form  is: 


^rs  ~ 


3  3  3  6  3  6  ^i2  1 

a=li,==ln=l;;=lfi7=l;;y=l/=hl=li2=l  it=l/l=l/2=l  12  “’2/ 


(5.18) 


where  the  powers  of  the  x  and  y  terms  in  the  area  integral  are 

.b,s 


t-  /•!  f-T 


"■«  =  "'^iUJ  *"1  *"‘n  *"'12 

"rs  =  "ftQi  *"4i\jjj  *4  *"n  *  4i 


(5.19) 

All  elements  of  the  stiffness  matrix  are,  thus,  linear  combinations  of  the  same  family  of  integrals 
over  the  panel's  area  of  the  form 


hR  (ra,n)  “  >’ 


(5.20) 


Note  the  explicit  dependence  on  thickness  coefficients  and  fiber  directions.  The  shape 
dependence  is  more  complex.  The  coefficients  Ui  and  Vj  (Eq.  3.7)  depend  on  the  x,y  positions  of  the 
vertices  of  the  panel.  These  coefficients,  in  turn,  determine  the  F  coefficients  in  Table  5.1.  In  addition, 
the  area  integrals  (Eq.  5.20)  depend  on  the  shape  of  the  panel  through  the  limits  of  integration. 


CHAPTER  6 


GEOMETRIC  STIFFNESS  MATRIX 

6. 1  Introduction 

Equations  for  the  elements  of  the  geometric  stiffness  matrix  are  derived  in  this  chapter  based  on 
in-plane  loads  from  a  wing  box  structural  analysis.  It  is  shown  that  the  geometric  stiffness  matrix  can  be 
expressed  as  combinations  of  area  integrals  of  simple  polynomials. 

6.2  Wing  Box  Stress  Analysis  and  In-nlane  Loads 
The  geometric  stiffness  matrix  was  derived  in  Chapter  2  as: 

(6.1) 

area 


The  [Kg]  matrix  thus  depends  on  the  derivatives  of  the  admissible  functions  contained  in  [F2],  and  the  in¬ 
plane  loads  contained  in  [N].  In  classical  linear  analysis  of  panels  the  in-plane  loads  are  assumed  given. 

In  the  case  of  wing  structures,  in-plane  loads  can  be  based  on  wing  box  stress  analysis  from  standard  finite 
element  techniques.  In  this  case,  in-plane  forces  along  the  sides  of  the  panel  are  obtained  by  some 
functional  approximation  based  on  either  nodal  forces  acting  on  nodes  surrounding  the  panel,  or  in-plane 
stresses  evaluated  for  the  finite  elements  surrounding  the  panel.  Alternatively,  based  on  stress  distribution 
in  the  panel,  in-plane  loads  for  buckling  analysis  can  be  evaluated  throughout  the  panel,  to  be  integrated 
over  the  area  of  the  panel  to  obtain  the  matrix  [KgJ.  (Refs.  43-50), 

For  preliminary  design  purposes,  if  the  skin  panels  are  small  relative  to  the  wing,  flutter  or 
buckling  analysis  may  be  accurate  enough  if  average  Nx,  Ny  and  Nxy  are  used  for  the  panel.  These  in¬ 
plane  stresses  are  assumed  constant  throughout  the  panel.  This  simplifies  the  integrations  in  Eq.  6.1,  and 
makes  it  possible  to  use  interaction  formulas  for  fast  approximate  buckling  analysis  (Refs,  45,47  and  49) 
and  33). 

When  an  equivalent  plate  modeling  approach  is  used  for  the  wing  box  analysis  (Refs.  32-34),  the 
in-plane  skin  stresses  are  obtained  from  the  wing  generalized  displacements  calculated  in  the  wing  box 
stress  analysis  stage.  In  the  formulation  used  in  Refs.  32-34  admissible  Ritz  functions  for  the  wing  box 
analysis  are  given  as  polynomials  in  x  and  y  ,  the  global  coordinates  used  to  define  the  planform  of  the 

wing  and  the  shape  of  all  panels  (Fig.  3.2).  In  this  case,  the  transverse  displacement  of  the  wing  is 
_  fsf-W  — .  — _ 
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where  a  bar  associates  variables  with  the  wing  box  analysis.  The  powers  m,  ,  n,  ,  and  the  number  of 

terms  N  w  are  known  from  the  Ritz  series  used  for  the  wing  box  displacement  solution.  The  coefficients 
q  i ,  are  the  generalized  wing  box  displacements.  Multiple  load  cases  can  be  accounted  for  in  the  wing 
analysis,  leading  to  different  generalized  displacement  vectors  {q} . 

In  equivalent  plate  wing  structural  analysis  based  on  classical  plate  theory  using  Kirchoffs 
kinematics  for  a  wing  with  a  symmetric  cross  section,  the  engineering  strains  in  the  x  and  y  directions  are 
given  by  (Refs  32-34) 


—  _  3  V  _  — 

^-“zW.yy  (6.3) 


Yxy  ^^zWjcy 


Let  the  wing  depth  be  given  by  H(^x,y) .  Then,  when  skins  are  thin  compared  to  the  depth,  they 

H  / 

can  be  assumed  located  at  z  =  ±  on  the  upper  skin,  in-plane  strains  are: 


Ex 

'w,xx  ' 

_  H(x,y) 

2  ' 

'W.yy-^ 

lW,xy^ 

Now,  in  a  polynomial  based  formulation  for  the  wing  box,  the  depth  of  the  wing  is  given  in 
polynomial  form 

—  —  {fnh-u) 

H{x,y)  =  I  Hih  x  ih  .y^  ih'  (6.5) 

ih=\ 

Since  the  displacement,  W  (x,y)  and  the  depth,  H(x.y) ,  are  polynomial,  it  is  evident  (Eq.  6.4), 
that  skin  strains  due  to  wing  deformation  are  polynomial  too.  There  are  total  of  Nl  layers,  each  with  fiber 
direction  0, ,  and  thickness  as  described  by  a  polynomial  equation.  The  in-plane  skin  stresses  in  each 
layer  of  the  panel  are  obtained  from  the  strains  by  the  constitutive  law: 
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where 


peep] 


was  derived  in  Chapter  5  as: 


p(6/)]=  [[2o]  +  [ei]cos2e^-  +[e2]cos4e  •  +[e3]sin2ej-  +[e4]sin4ej] 

The  in  plane  loads  can  now  be  found  in  terms  of  wing  box  displacements  by  integration  through 
the  thickness  of  the  panel: 


A^J=j  a,U=J[e,e)^.  U=- 


rf^ 


2(9)  Jfe 


The  in-plane  stiffness  matrix  is  defined  as 

[^]  =  J  [2(z)]^fe 

z 

which  was  derived  in  Chapter  4  to  be 

r=  ^  inii )  (ni ) 

[A]=  .s  I^[e(e,)]  x  k  .y  k\Ti 

which  leads  to 


In  order  to  express  Eq.  6.11  in  terms  of  the  generalized  displacements  {q}  calculated  in  the 
wing  box  solution,  we  define  a  new  vector  {OC }  of  wing  box  curvatures: 

\W^\  ...If  ■  ■ 

{a}='  Vr,yy  *=  ...  npw^pw  -  ^  ... '-  •  i6.l1 

2W  XV  ~  ~  mpw-l  npw-l  ^ 

,xy  j  mpwnpwx  y  ^ 
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where  pw  is  an  index  for  the  p’th  element  in  the  wing  Ritz  series.  There  are  N  w  terms  in  the  Ritz  series 
for  wing  displacement.  The  matrix  [a]  is  thus  3  by  N  w.  Each  element  of  this  matrix  is  of  the  form 

^qWypw  '  ^ ^  the  coefficient  Wq^  py^  ,  and  powers  of  x  and  y, 

^qw,pw  ^qw,pw  respectively,  are  described  in  Table  6.1. 

Table  6.1  Coefficients  and  Powers  of  x  and  y  of  the  Matrix  [a] 


^qw,pw 

'^qw,pw 

^qw,pw 

qw=l  (Row  1  of  the 

m  pw(m  pw’l) 

m  pw-2 

npw 

Matrix) 

qw=2  (Row  2  of  the 

npw(npw-l) 

m  pw 

npw'’2 

Matrix) 

qw=3  (Row  3  of  the 

m  pwH  pw 

m  pw-1 

npw-1 

Matrix) 

When  the  index  qw  denotes  the  row  of  the  matrix  (its  values  varying  from  1  to  3)  and  the  index 
pw  specifies  terms  of  the  polynomial  displacement  series  for  the  wing  box  (varying  from  1  to  N  w),  then 
the  in-plane  load  Nx  (Eq.  6.11)  for  the  panel  can  be  expressed  in  polynomial  form  as 


H  3 


H  3 


Nw 


2  o  ^  ^  ,'^qw,pw 

I  qw=\  ^  2.  cfw=\  pw=l 


'^\jw,pw  ^qw,pw 


^  pw 


leading  to 

Ji  3  Nw 


J^,pw  P'qw,pw- 

y  H  pw 


(6.13) 


Similarly  the  expressions  for  Ny  and  Nxy  can  be  derived; 


H  3  Nw 


/\w,pw  J^qw,pw 


^  pw 


(6.14) 
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H 


3  Nw 


Nxy=---2  X  X,A, 


2  V=i;^=i  ^ 


’'^\]w,pw  ^qw,pw  ~ 


^  pw 


—  3  Nw 

“V=lj=/34>»‘V.PX' 


(6.15) 


J^,pw  J^qw,pw- 

^  y  H  pw 


Now  the  polynomial  expression  for  wing  depth  (Eq.  6.5)  can  be  substituted  into  Eqs.  6.13-6.15. 

r  .  N. 

The  general  expression  for  terms  of  the  [N]  matrix  ,  [A(  J  = 


xy 


N„.  N.. 


,1S 


N 


ppm 


I  3  Nw 


(mif^+mqw,pw)  (nif^+nqw,pw)- 

y  ^  pw 


(6.16) 

The  indices  pp  and  qq  can  be  either  1  or  2.  Note 
If  pp  =  1  and  qq  =  1,  then  Appp 

If  pp  =  2  and  qq  =  2,  then  A  ppp  =  A  2  q^  (6.17) 

Whenpp^fqq,  Appp  q^  =2A3  q^ 

The  polynomial  expression  for  the  [A]  matrix  (Eq.  6.10)  is  now  used: 


^pp^q  - 


2^  3 

2ih^=lp^=lSlkll^^PPP’^ 


■^ih 


^  pw 


(6.18) 


where  the  index  ppp  means  the  following: 

If  pp  =  1  and  qq  =  1,  then  Qppp^qw  “  Ql,qw 

If  pp  =  2  and  qq  =  2,  then  Qppp,qw  =  Q2,qw  (6.19) 

When  pp  ^  qq,  Qppp^qw  =  ^Q3,qw 

The  matrix  [N],  then,  can  be  expressed  as  a  polynomial  in  x  and  y  according  to  Eq.  6.18.  This 
equation  shows  how  [N]  depends  on  the  wing  box  solution  ,  the  depth  of  the  wing,  the  thickness 
coefficients  for  layers  in  the  panel,  material  properties  and  fiber  directions.  Of  course,  the  wing  solution 
{q} ,  also  depends  on  depth,  thickness  of  layers,  fiber  directions  and  material  properties.  These  affect  the 


j 


stiffness  matrix  of  the  wing  as  described  in  Refs.  32-35.  It  should  be  emphasized  again  that  polynomial 
in-plane  loads  for  the  flutter  analysis  can  be  obtained  in  a  similar  manner  from  wing  box  analysis  based 
on  first  order  shear  deformation  plate  theory  (Refs.  35  and  37)  or  from  finite  element  results,  when  skin 
stresses  are  approximated  by  polynomials  using  least  square  fitting  (Ref.  39) 


6.3  The  IF;!  Matrix 

The  [Fa]  matrix  was  defined  in  Chapter  2  as 


/l.y 


fl.x 


f N,x 

f N,y 


(6.20) 


The  pth  admissible  function  is  defined  in  Chapter  3  as: 


3  6 


The  first  derivatives  are  simply 


=1; 

3  6 


(m^+mp  (n^+n^+n}--l) 


■y 


(6.21) 


(6.22) 

(6.23) 


Thus,  the  element  q,p  of  the  matrix  [F2]  is  of  the  form: 


3  6  ...  (mflfp  {nflfP) 

Fo  =  I  I  F^9’P-x  -y 

2(q,p)  ^i,J  ^ 


(6.24) 


The  coefficients  and  powers  of  x  and  y,  mflf'f  and  ,  defining  the  elements  of  the 

[F2]  matrix,  are  described  in  the  table  6*2. 
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Table  6.2  Coefficients  and  Powers  of  Terms  of  the  [F2]  Matrix 


p  ^■>p 
^2iJ 

mflff 

nnfj 

q=l  (Row  1  of  [F2] 

UiVj(mp+mJ) 

m  p  +  iMj  -1 

n^  +  n“+nj 

Matrix) 

q=2  (Row  2  of  [F2] 

UiVj(n;'  +  n“+n[) 

W  V 

nip +mj 

W  U  V 

Dp  +ni  +nj -1 

Matrix) 

6.4  The  Geometric  Stiffness  Matrix 

Elements  of  the  matrix  [Kq]  can  now  be  expressed  in  polynomial  form  as  follows: 


2  2 

^  a=lb=l 


^2  a.,r^ajb^2 


b,s 


dxdy 


(6.25) 


The  indices  r  and  s  identify  terms  in  the  panel  Ritz  displacement  series  for  the  panel.  Using 
polynomial  expressions  for  [F2]  and  [N]  (Eqs.  6.18  and  6.22),  the  r,s  element  of  [Kq]  is  written  as: 


N 


K. 


Nl  n,. 


Gr,& 


I  2  2  3  6  3  6  3  Nw 

^  a=lb=li=l  j=lii=l  jj=lih=lqw=l pw=l  it=lk=l 


p  a,r 


p  b,s 


(6.26) 


■Qppp,qw  (e.p  •  ^qw,pw  ■ 


■t: 


it  ^  cc  (nGp  g)  j  j 

■<lpw]}X  y  dxdy 


where 


^^r,s  ^  +K‘ 


,  =nf2‘l''  +nf2\:^..  +%  +nf'‘ 


r,s 


The  index  ppp  used  in  Q  pppqw  (0;^ )  defined  in  Eq.  6.19.  As  in  the  case  of  the  mass  and 

stiffness  matrices,  the  geometric  stiffness  matrix  is  represented  as  summation  of  surface  integrals  of 
polynomial  terms  calculated  over  the  area  of  the  panel.  Integrals  of  the  same  family 


^TR  (m,n) 


Wx^y^dxdy 

Area 


(6.27) 
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are  used.  Evaluation  of  these  integrals  is  accomplished  using  die  same  subroutines  used  in  the  wing 
analysis  and  all  panel  matrix  analysis. 

6.5  Further  Discussion  of  Integration  With  Wing  Structural  Analysis 
Recall  that  the  results  of  the  wing  analysis  affect  analysis  of  the  panel  through  the  in-plane  load 
matrix  [N].  This  matrix  is  obtained  by  integrating  skin  stresses  (due  to  wing  deformation)  through  the 
thickness  of  the  skin,  and  the  resulting  Nx,  Ny  and  Nxy  are  functions  of  x  and  y  as  evident  in  Eqs.  6.18 

and  6.19.  Now,  while  equivalent  plate  wing  structural  analysis  leads  to  good  skin  stress  predictions  in 
stiff,  low  aspect  ratio  wings  (Refs.  32-34),  still  equilibrium  is  not  guaranteed  for  isolated  skin  panels. 

Thus,  the  in-plane  loads  (as  given  by  Eq.  6.18)  are  generally  not  in  equilibrium  (This  problem  is  not 
encountered  when  wing  analysis  is  based  on  the  finite  element  method,  when  in-plane  loading  on  the 
circumference  of  an  isolated  panel  is  determined  from  nodal  forces  on  that  circumference). 

One  way  to  avoid  this  difficulty  and  simplify  the  geometric  stiffness  matrix  is  to  use  an  average 
in-plane  loading  to  be  assumed  constant  over  the  panel  (Ref.  49).  If  point  (xq,  yo)  inside  the  panel  is  used 

to  evaluate  these  average  in-plane  loads,  then  the  elements  of  [N]  become 

N  — 

3  W 

^ pp,qq  ~ ~2’ ^_^^PPP^q^^q^ ^p^^q pw 

(6.28) 


The  expression  for  elements  of  the  geometric  stiffness  matrix  (Eq.  6.24)  has  to  be  modified,  since 
terms  including  xq  and  yo  become  constant  for  the  area  integration,  and  can  be  taken  out  of  the  integral. 


^  1223636  3  Nw  ^  ~  ^  ^ 

Qppp,qw  (6 .j. )  ■^qw,pw  ‘  ‘  '^jl  <lpw ' 

{mf  2?’r +m/2^.’'^.  .+mPJ)  +nf2^.'^.  ■+nP^) 

jjx  ^  -y  ^  dxdy 


(6.29) 


CHAPTER? 


AERODYNAMIC  FORCE  MATRICES 
7.1  Introduction 

The  explicit  expressions  for  the  terms  of  the  aerodynamic  force  matrices  are  derived  in  this 
chapter.  It  is  shown  that  the  aerodynamic  damping  and  stiffness  matrices  can  be  formulated  as 
combinations  of  area  integrals  of  simple  polynomials.  First  order  linear  piston  theory  aerodynamics  is 
used  for  the  formulation  of  the  matrices.  At  the  conclusion  of  the  chapter  the  characteristics  of  the 
aerodynamic  matrices  are  discussed. 


7.2  Piston  Theory 

As  shown  in  Chapter  2,  the  aerodynamic  contribution  to  the  panel  aeroelastic  system  is  given  by 


(7.1) 


where  q(x,y)  is  the  change  in  pressure,  Ap.  First  order  linear  piston  theory  (Ref.  1)  approximates  Ap  by 


(7.2) 


where  p„,  U«,  and  M„  are  the  free  stream  density,  velocity,  and  Mach  number,  and  ^  is  the  flow  direction. 


U 


\ 


^  Ml 


Figure  7.1  Panel  in  Angled  Flow 


The  simple  relation  between  directional  derivatives  in  the  (^,11)  and  (x,y)  planes  is 
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—  =  C0S(t)  — +  sin(j)-:— 
dq  dx  dy 


Simplifying  Eq.  7.2,  piston  theory  can  be  rewritten  as: 


Ap  =  cos(t)W;  +  sin(t)lV‘  +P,Wl 


where 


(7.3) 


(7.4) 


(7.5) 


^  -  1  f/ 


(7.6) 


7.3  Aerodynamic  Force  Matrices 

To  express  the  aerodynamic  forces  in  terms  of  the  generalized  coordinates  the  panel  displacement 
terms  and  their  derivatives  are  expressed  as: 


^  ~  ■|/l{Ar,y)  fl(x.y) 


W^xU>y)=|/l,x  /2,x  -  /w,x| 


N{x,y) 

^2 

Qn 


=[f.M 


'rj(x,y)=j/,,,  ... 


(7.7) 
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V/M,y)  =  \f,  f,  ...  /,|< 

Substituting  Eqs.  7.4  and  7.7  into  Eq.  7.1  gives: 

{e}=  ^^C0S(l)jL[F,f  ^^y  +  /’^sin(|)Jj  [Fif  dxdy  + 

‘>^y  <7  8) 


^N,t 


=  ['7,  ]{«.,} 


It  is  natural  at  this  point  to  break  up  the  aerodynamic  force  into  an  aerodynamic  stiffness  matrix 
and  an  aerodynamic  damping  matrix. 

{Q}  =  [QsHff  }+  Pt  [Qdaynp  }  (7-9) 

where 


Qst^  rs 


JJarea 


dxdy 


(7.10) 


Qdamp  =  Jf  Irfs^^dy 

^  dJ  area 


(7.11) 


fr  and  fs  are  the  r^  and  s^  admissible  functions.  In  chapter  3  the  admissible  functions  were  derived  in 
terms  of  the  coefficients  Ui  and  Vj  as 

3  6 


Pix,y) 


(7.12) 


i=l  i-1 


The  first  derivatives  of  the  admissible  functions  with  respect  to  x  and  y  are 

4..  =  i  i  K + 

i=l  J=i 


(7.13) 


39 


p.y 


i=l  ;=1 


(7.14) 


Substituting  these  into  Eqs.  7.10  and  7.11  shows  the  final  expression  for  the  elements  of  the 
aerodynamic  damping  and  stiffness  matrices. 


3  6  3  6 


1  =  1  ;  =  1  »  =  1  jjzzl 

+n-j  +n:)ll^^x-y-'dxdy 

i=l  y=l  n=l  jrysl 

e._,  „  =tttt^yiV,V,jl^_^x-y’dxdy 


/=!  y=l  ii={  jj=i 


where 

m  =  (mj  +  nCjj  +  rri^  +  ) 

n  =  «+<+n;+n; +«;+«;) 


(7.15) 

(7.16) 


The  elements  of  the  aerodynamic  stiffness  and  damping  matrices  are  linear  combinations  of  the 
same  integrals  as  the  mass,  stiffiness,  and  geometric  stiffness  matrices. 

The  aerodynamic  stiffness  matrix  is  dependent  on  the  air  stream  angle,  (j),  the  dynamic  pressure 
and  Mach  number  through  P^,  and  the  panel  shape  variables  through  Uj  and  Vj .  The  aerodynamic 
stiffness  matrix  is  dependent  on  the  dynamic  pressure  and  Mach  number  through  Pt,  and  the  panel  shape 
variables  through  Uj  and  Vj.  Of  course  the  area  integrals  themselves  depend  on  the  panel’s  planform 
shape  variables. 

In  the  stability  analysis,  as  shown  in  Chapter  2,  it  is  convenient  for  the  eigenvalue  problem  to 
express  the  aerodynamic  matrices  in  the  form: 


U  1= _ 


(7.17) 


which  leads  to  an  expression  for  the  aerodynamic  generalized  forces  in  the  form: 
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{Q)=P.V-[A^\q,}+P.Ul [ V  ]{^ } 


(7.18) 


7.4  Characteristics  of  the  Aerodynamic  Matrices 

Both  Qstiff  and  Qdamp  have  interesting  characteristics.  Looking  at  Eq.  7.11  it  is  obvious  that  Qdamp 
is  a  symmetric  matrix  -  a  fact  which  allows  faster  computation.  As  a  matter  of  fact,  if  the  panel  is  made  of 
a  single  material,  and  has  a  constant  thickness,  the  aerodynamic  damping  matrix  is  proportional  to  the 
mass  matrix. 

As  to  the  matrix  Qsuff  -  with  the  current  choice  of  admissible  functions  (for  the  simply  supported 
case)  this  matrix  is  a  skewed  synunetric.  In  addition,  certain  elements  in  this  matrix  are  identically  zero. 

The  skew-symmetry  of  Qsuff  is  shown  by  applying  integration  by  parts  with  respect  to  x  and  y  to 
Eq.  7.10: 


Qstiff  rs  =C0S(j)j  -jfj.jx  cfy-t-sin(t)J 


frfsC-jfsfr.ydydx  (7.19) 


Because  the  admissible  functions  are  zero  along  the  boundaries,  the  boundary  terms  are  zero  leaving 


Qstff  n  dxdy 


Clearly: 
Qstiff  rs 


=  -Q 


stiff  sr 


(7.20) 


It  can  also  be  shown  that  several  terms  including  all  the  diagonal  terms  in  the  aerodynamic 
stiffness  matrix  are  zero.  Starting  from  Eq.  7.10  with  r=s: 


a,  „  =cosc|iJL/./,,  dxdy  +  sin(^llj,ff, 


dxdy 


(7.21) 


Looking  at  the  integrands: 


and 
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2d  y 


(7.22) 


Therefore  any  diagonal  term  is  identically  zero.  Of  course,  the  diagonal  elements  of  any  skew-symmetric 
matrix  are  zero.  However,  the  derivation  above  serves  to  illustrate  how  other  terms  in  the  present 
formulation  may  be  identically  zero.  For  example,  if  the  flow  is  directed  along  the  x  axis,  sin(t)=0,  two 
admissible  functions  may  be 

and 

fs=fB(x,y)-x'”y'‘ 

Thus,  in  effect  for  this  special  case 

/.  =  /.•/ 

Now,  the  area  integral  associated  with  the  r,s  element  of  the  aerodynamic  stiffness  matrix  is 
JJ  fr  ■  fs^dxdy  =  jj/^  •  ■  y^cbcdy  =  ^  y^dxdy 

Integration  in  the  x  direction  will  immediately  show  that  in  the  case  of  simply  supported  panels  (fB(x,y)=0 
on  the  boundary)  this  area  integral  is  zero.  Thus,  using  complete  polynomials  as  the  functions  to  multiply 
fB(x,y)  some  fr  and  fs  have  the  same  powers  of  x  or  y  for  r  For  these  terms,  and  depending  on  the  flow 
direction,  the  x  or  y  derivative  components  of  Qsmfrs  will  also  be  zero. 


CHAPTERS 


ANALYTIC  SENSITIVITIES 

ST  Introduction 

Expressions  for  the  analytical  sensitivities  of  the  system  matrices  with  respect  to  the  design 
variables  are  derived  in  this  chapter.  Eigenvalue  sensitivities  and  flutter  speed  sensitivities  are  also 
developed, 

S.2  Mass  Matrix  Sensitivities 

The  mass  matrix  is  dependent  on  the  panel  thickness  design  variables  and  the  shape  design 
variables.  The  explicit  expression  for  the  mass  matrix  (Eq.  4.9)  allows  for  easy  determination  of  these 
analytical  sensitivities. 


The  mass  matrix  is  dependent  on  Tj[ ,  the  sizing  variable  corresponding  to  the 
layer,  through  the  vector  T,^  (Eq.  4.9).  The  derivative  of  a  mass  matrix  term  is  given  as: 

3  d  a  f, 
a  Tl  dfkdTl 

Recall  T|^  is  defined  as: 

i=l 

Therefore,  for  a  given  layer,  i,  and  thickness  term,  k, 

d  fk  _ 

Substituting  this  back  into  Eq.  8.1  shows  that 

3  T‘  3  T, 

for  every  layer  i.  Differentiating  Eq.  4.9  for  any  specific  k  and  i  gives: 

a  A/„  a  3  6  3  6 


where 
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m  =  m]  +  m;;  +  nC  +  +m[ 

J  JJ  r  s  K 

n  =  n-  +n,“  +n^.  +<  +n^  +n' 

Notice  that  if  all - = —  are  calculated  before  the  mass  matrix,  then  the  mass  matrix  can  be 

d  Tk 


formed  as: 


This  allows  for  the  mass  matrix  and  its  thickness  sensitivities  to  be  calculated  at  the  same  time  with  no 
additional  computations  that  would  be  required  by  calculating  the  mass  matrix  alone. 


8.2.2  Mass  Sensitivities  With  Respect  to  Shape  Design  Variables  Vt  .  vg.  x^.  and  x^ 

The  mass  matrix  is  dependent  on  the  shape  variables  through  the  coefficients  Ui  and  Vj  (Eq.  3.7) 
and  the  integrals  lTR(m,n)  because  the  limits  of  integration  depend  on  the  shape  variables  (Eq.  4.10), 
Analytical  sensitivities  for  Ui  and  Vj  are  obtained  through  direct  differentiation  of  Eq.  3.7: 

Derivatives  with  respect  to  yt 

d  u, 

^  Yl 

a  u. 


=  yF 


=  -i 


=  0 


^  Yl 

a  U3 
a  yl 

a  Vl  ^AlYr  ~XarYr  p  ■  ^flYr  ~  ^frYr  o 

(y.-yj  '  (y.-y0“  ‘ 

a  V2  _  ^AlYr~^ArYr  ^FlYr~XfrYr 
^  (Yr-Yl)'  (Yr-Yl)' 


(8.7) 
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^  ^3  _  ^FlYr  ^FrYr  c  _l  ^AR  ^AL  d  .  ^ALyR  ^ArXrc  ,  ^FR“^FL  r. 

^  yi-  {y^-yx)  (yr-Yl)  (yr-Yl)  (yr-Yl) 


a  V, 


^  Yl 

a  V, 


4  _ 


=  0 


^fr  ^FL 


^AR  ^AL 


^y^  (Yr-Yl)'  (Yr-Yl)' 

^  ^6  _  ^FR  ~  ^FL  o  I  ^AR  ~  ^AL  c 

3y.'(y.-yJ  ‘  (y.-y.)'  ' 


Derivatives  with  respect  to  Yr 

a  u, 

a  Yr 


a  u, 

5  Yr 

a  U3 

^  Yr 

^  ^1  _  ^ArYl  ~^AlYl  r»  ,  ^FrYl  ~  ^FlYl  r> 

^Yr"  (y,-yj  ^ 


=  -l 


=  0 


(Yr-Yl) 

^  V2  ^  x^RyL-XALyL  Xf^y^-XpLyi, 

^  (yr-yl)' 


(Yr-Yl) 


(8.8) 


^  ^3  _  ^FrYl  ^FlYl  o  ,  ^AL  ^AR  n  ,  ^ArYl“^AlYlc  ■  ^FL“^FR  r. 

y2  ^A  +  ,  x2*^F+  ,  .,2  ^P+-  :j^A 


^y^  (Yr-Yl)' 


(Yr-Yl) 


(Yr  -  Yl) 


(Yr-Yl) 


a  y 


^  Yr 

a  V 


4  _ 


=  0 


ypL  ^FR 


Xal-Xar 


^  yR  (Yr-Yl)'  (Yr-Yl)' 

^  '^6  _  ^FL  ~  ^FR  c  .  ^AL  ~  ^AR  c 

^yR'(YR-YL)^^  (Yr-Yl)^^ 
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Derivatives  with  respect  to  Xal 

Mi_  =  o 

3  U2 
d  Xal 

d  Xal 


=  0 


a  V, 


Yr 


a  X 


AL 


a  V, 


Yr-Yl 

Yr 


R. 


a  X 


AL 


a  V, 


Yr  -Yl 
-1 


Rp  +- 


Yr 


^  Xal  Yr-Yl  Yr-Yl 


a  V, 


a  X 


=  0 


AL 


a  v< 


5  _ 


1 


a  X 


AL 


a  V, 


Yr  -  Yl 
-1 


a  X 


AL  Yr-Yl 


Derivatives  with  respect  to  xar 


a  Ui 

=0 

d  Xar 

a  u. 

=0 

1 

^  Xar 

a  U3 

=0 

a  Xar 

a  % 

_  -Yl 

AR  Yr-Yl 


a  X 

A3.-  yi 


R. 


(8.9) 


(8.10) 
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a  V. 


1 


a  X 


AR 


Yr  -Yl 


•Rp  +■ 


-Yl 


Yr  -  Yl 


a  y 


4  _ 


a  X 


=  0 


AR 


a  V, 


5  _ 


-1 


a  X 


AR 


Yr  -  Yl 

1 


a  V. 


9  Xar  Yr  -  Yl 


Derivatives  with  respect  to  Xfl 

a  u, 


a  X 


=  0 


FL 


a  u 


^  _ 


a  X 


=  0 


FL 


a  u 


3  — 


a  X 


=  0 


FL 


3  V|  Vr 


a  X 


FL 


a 

a  Xr 


Yr-Yl 

-Yr 


^FL  Yr-Yl 

^  ^3  ^  Yr 
^  XpL  Yr-Yl 

a  V 


Sa  +■ 


-1 


Yr-Yl 


4  _ 


a  X 


0 


FL 


a  V 


5  _ 


1 


a  X 


FL 


a  V. 


Yr-Yl 

-1 


a  X 


FL  Yr-Yl 


(8.11) 
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Derivatives  with  respect  to  Xfr 

a  u, 


a  X 


=  0 


FR 


a  u. 


a  X 


a  u. 


a  X 


a  V, 


a  X 


a  V, 


a  X 


FR 


a  V. 


a  X 


a  V, 


FR 


a  X 


a  V, 


a  X 


0 


FR 


a  Xpj^ 
a  V, 


=  0 


R. 


FR 


FR 


Yr-Yl 
Yl 

Yr-Yl 

_  -Yl 
Yr-Yl 

=  0 


Sa  +■ 


1 


R. 


Yr-Yl 


FR 


5  _ 


1 


FR 


Yr  -  Yl 

Yr-Yl 


(8.12) 


r™  ^  ^TR  _ 

The  derivatives  ,  where  x  is  any  shape  variable,  are  prepared  using  Refs.  34  and  35. 

O  X 

The  shape  sensitivities  of  the  area  integrals  are  linear  combinations  of  other  members  of  the  same  table  of 
integrals.  No  new  integrations  are  needed  for  the  sensitivities. 

With  this  information,  the  derivative  of  the  mass  term  Mrs  is  calculated  as  follows: 


dM 

d  X 


6  3 


k=\ 


i=l  ;=1  «=1  jj=\ 


\d{u,V,U.Vj,), 

'  dx 


TR(m,n) 


d  X 


with 

m  =  mj  +  m'^j  +  +  mj  +  ml 

n  =  n“  +n,"  +n'’j  +n‘).  +<  +nl 


(8.13) 
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8.3  Stiffness  Matrix  Sensitivities 

With  tlie  explicit  expression  of  stiffness  matrix  elements  in  terms  of  thickness,  fiber  directions 
and  shape  of  the  panel  available  in  Eq.  5.18,  it  is  straightforward  to  obtain  sensitivities  analytically. 


8.3.1  Stiffness  Sensitivities  With  Respect  to  Thickness  Design  Variables  t|- 

Planform  shape  variables  and  orientation  angles  are  fixed  in  this  case.  The  thickness  coefficients 
T  appear  in  the  expression  for  the  stiffness  matrix  (Eq.  5.18)  explicitly  in  a  triple  summation  over  the 
indices  k,  11  and  12.  Sensitivity  is  then  obtained  by  direct  differentiation,  noting  that  if  the  design  variable 

Q  ar*  / 

involved  is  ,  then  =  1  only  when  i=q  and  k=r.  Otherwise,  the  derivative  is  zero. 

/ar^ 

/  r 

Differentiating  Eq.  5.18  leads  to: 

BK  3  3  3  6  3  6  ^i2  I  nr  -  /  \ 

- IIIXEIXSI  III  —P,VnPi/inQaJb{^i\ 

3  7^  a=lfc=l/i=l  jf7=liH=l ;>>■=!  i=in=l/2=l  it=l/l=l/2^1 12  ^*>2/  «/ 


"‘rs  =  "'fufjj  *4  *"'tl  *”l2 

"rs  ""ftifij  *4  *4\  +";2 


8.3.2  Stiffness  Sensitivities  With  Respect  to  Fiber  Direction  9 . 

The  angle  represents  the  direction  of  fibers  in  the  i-th  layer,  and  the  stiffness  matrix  depends 

on  0 .  through  elements  of  the  matrix  [ece  ■ ) I  (Eq.  5.6).  The  derivative  of  each  term  in  the  stiffness 
matrix  will  be  calculated  in  the  following  manner.  In  the  summation  (Eq.  5.18)  over  i=l  to  Nl,  all 
matrices  ^(6 , )  j  are  set  to  zero,  except  the  matrix  ^(0 , )  j  corresponding  to  the  0 variable 

considered.  This  particular  |c(0 , )  j  is  replaced  (in  Eq.  5.18)  by  the  following  expression: 
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=  -2[(2i]cos26^-  -4[g2]cos40^-  +2[(22]sin20j  +4[(24]sin46j-  (8.15) 

i 

Equation  5.18  is  ,  thus,  used  for  the  sensitivity  of  the  stiffness  term,  with  the  derivative  Eq.  8.15  replacing 


8.3.3  Stiffness  Sensitivities  With  Respect  to  Panel  Planform  Variables  vt_.  vp.  xpp,  x^i  and  xar 
Thickness  coefficients  and  orientation  angles  are  held  fixed.  The  terms  Krs  of  the  stiffness 
matrix  depend  on  the  shape  through  the  matrix  [F3]  and  the  integrals  IiRdn^i)  •  If  x  is  any  planform  design 
variable  then: 


d  u.  a  Vj 

- -V:+U,  ' 


'  a 


V 


(Integers  from  table  5.1) 


(8.16) 


The  derivatives  of  the  Ui  and  Vj  coefficients  and  the  integrals  lTR(in.D)  are  the  same  as  for  the  mass  matrix 
developed  in  section  8.2.  With  this  information  the  derivatives  of  the  stiffness  matrix  can  be  calculated 
as: 


a  Kfs  _  3 
ax 


3  3  6  3  6  ^il  ^i2  l 

YYY  — 

a=lZ7=ln=l  jj=liihl i/^1 /=lil=lz2=l  it=l/lti/2=l  12 


ar.^’C. 

II  ^jj 


dFb^... 

pu,i  ,  pCl,V  m^JJJ 

iiijjj  iijj' 


nb,s 


d\ 


Qa^(^i yi  ■  r//  •  TI2  ■  Ijp  (8. 17) 


,ra,r 


9x 


8.4  Geometric  Stiffness  Matrix  Sensitivities 

Equation  6.24  gives  the  geometric  stiffness  matrix  explicitly  in  terms  of  the  thickness,  fiber 
direction,  and  shape  design  variables. 

8.4.1  Geometric  Stiffness  Sensitivities  With  Respect  to  Thickness  Design  Variables  Tj^ 

The  design  variable  in  this  case  is  the  k-th  coefficient  in  the  polynomial  thickness  series  for  the  i- 
th  layer.  Examination  of  Eq.  6.24  reveals  that  the  geometric  stiffness  matrix  is  explicitly  linear  in  the 
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thickness  coefficients  T  .  It  is,  of  course,  also  dependent  on  those  coefficients  via  the  wing  box  solution 
{q}  ,  unless  it  is  assumed  that  in-plane  loads  [N]  do  not  change.  Differentiation  of  Eq.  6.24,  using  Eq. 


33  leads  to 


^%r,s  1  2  2  3  6  3  6  3  Nw  ^  ar  ~  bs 

ar. 


It 


Z  S  Z  X  Z  S  Z  Z  Z  'Z,  1,  F  F  . 

2a=lb=\i=lj=Ui=ljj=lih=lqw=lpw=l  R=IS=1  ^HJj 


k 

where 
^^r,s  " 

^^r,s  '^^ih  +"^Jt  +”yv.p». 

it  R 

and  is  equal  to  1  only  when  it=R  and  k=S.  Otherwise,  it  is  zero. 


(8.18) 


8.4.2  Geometric  Stiffness  Sensitivities  with  Respect  to  Fiber  Direction  6,^ 

Layer  orientations  affect  the  geometric  stiffness  matrix  through  the  material  matrices 
Qppp,qw  (0  .j.)  '’''“8  generalized  displacements  {q] .  The  analytic  sensitivity  with 

respect  to  fiber  direction  in  a  given  layer  is: 


^%rs  1  2  2  3  6  3  6  3  Nw  ~  ar~  h, 

^  a=\b=\i=\j=Ui=\  jj=\ih=\qw=\pw=\  itt=\.k=\  ^b,jj 

^Qppp^iQ:„)-  =  dqpyy  .. 

[  _  qpw  +Qppp^  ae. 


ae, 


It 


It 


Where  the  derivatives  of  Qppp^qw  (0  are  the  same  as  for  the  stiffness  matrix  (Eq.  8.15) 


(8.19) 
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8.4.3  Geometric  Stiffness  Sensitivities  With  Resnect  to  Planform  Variables  vi^.  vp.  xfi_.  xfr^  ai  and 
XAR 

Thickness  coefficients  and  orientation  angles  are  held  fixed.  The  geometric  stiffness  matrix  Kq 

depends  on  the  planform  variables  through  U,  and  V,  terms  in  F of’f  and  F (Table  6.2  and  Eq. 

6.24).  There  is  also  a  dependence  on  the  area  integrals  Ijr  as  discussed  previously.  The  derivatives  of 
the  [F2]  terms  are  calculated  by: 


d 

d  X 


(  \ 

at/.  a  v,. 

— i-v.  +t/. — 

d  X  ^  '  a  jc 


(Integers  from  table  6.2) 


(8.20) 


which  is  substituted  into  the  derivative  of  Ko, 
a  %rs  1  2  2  3  6  3 


6  3  Nw 


Jrs  o  /z  /vw  L 

j  ^90  _ _ ^  ^  ^ 

X  2  a=l  b=l  t=l  7=1  ii=l  jj=\  ih=\  qw=\  pw=\  it=\  k= 

a  . 

—b,s 


{( 


[q ppp,qw  (6  .p  •  ^qw,pw  '  f^ih  ’  ) '  | 


d  X 


2fi,;; 


+  F. 


'a,r 


3 


i,j  d  X 


'^pwhR(m,n) 


— a,r — b,s 


r 


'  piw 


^  V/?(m,n)  ^^pw 

d  X  ^  a  X  TR  (m,n) 


M 


(8.21) 

where  the  powers  of  integrands  in  lTR(m,n)  m=mGr,s  and  n=nGr,s.  In  Eq.  8.21  it  is  assumed  that 
overall  wing  planform  is  fixed,  and  panels  are  changing  shape  and  location  due  to  moving  of  control 
surfaces,  ribs  and  spars.  If  overall  planform  shape  of  the  wing  is  changing,  then  derivatives  of  the  wing 
depth  coefficients  with  respect  to  the  shape  design  variables,  d  jd  x  ,  must  be  added,  since  the  wing 

depth  is  defined  in  global  x,y  coordinates  (Eq.  6.5). 

When  a  constant  in-plane  load  matrix  is  assumed  using  loads  at  a  point  (xo,yo)  on  the  panel  (Eqs. 
6.27  and  6.28)  then  the  motion  of  point  (xo,yo)  must  be  taken  into  account  in  the  Kg  shape  sensitivity. 
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8.5  Aerodynamic  Matrix  Sensitivities 

The  aerodynamic  stiffness  and  damping  matrices  are  dependent  only  on  the  shape  design 
variables  through  Ui  and  Vj  and  the  integrals  lTR(mji).  The  sensitivities  are  obtained  by  direct 
differentiation  of  Eqs.  7,15  and  7.16. 


3a,» 


■stiff  rs  3  6  3  6  (U  V  U  V  ) 

1=1  ;=1  ,7=1  y7=i  [  ox 


d  X 


’jj  '  ^^TR{m-\,n) 

dl. 


+ 


‘  TR  (m  -l/t) 


d  X 


3  6  3 


sin(t)SSXS 

,=l  j={  ,7=1  jr)*=i 


d  X 


K+n;. +«:)/. 


TR  (m,n  -1) 


(V,V,U,V,iXn-  +n-j, 


d  X 


(8.22) 


damp  rs 


d  X 


,=1  j=l  ,7=1  yy  =  l  I  O  X 


TR  (m/i) 


^  y  .JJ  .y  .^^LbShiH. 
'  ^  dx 


(8.23) 


in  both  equations  m  =  (mj  +  m],  +  /n“  +  )  and  n  =  (n“  +  n,“  +  +  n  Jy  +  +  ^7 ) 

Recall  from  Eq.  7.17  that  the  coefficients  relating  Ajtiff  and  Adamp  to  Qsoff  and  Qdamp  are  functions 
of  the  Mach  number  only,  which  is  held  constant  in  this  analysis  (although  sensitivities  with  respect  to 
Mach  number  are  easy  to  calculate).  Therefore  the  sensitivities  of  Asaff  and  Apamp  are: 


dA 


■stiff 


aov 


dA 


^damp 


aov 


1  ^  Qstiff 

3DV 

1  Ml-2dA^f 

aov 


(8.24) 


8.6  More  on  the  Analytical  Sensitivities 

It  is  clear,  examining  the  analytic  sensitivities  of  the  system  matrices  derived  in  the  previous 
sections,  that  all  involve  area  integrals  of  polynomial  terms  over  the  panel’s  area.  Since  members  of  the 
family  of  integrals  lTR(m,n)  are  generated  for  the  analysis  stage,  there  is  no  need  to  generate  them  again 
for  the  sensitivity  calculation  stage.  No  numerical  integration  is  therefore  needed  for  either  analysis  or 
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sensitivity  calculations.  This  is  similar  to  the  process  by  which  analytic  sensitivities  are  obtained  for  the 
wing  box  analysis  (Refs.  34  and  35).  A  table  of  area  integrals  over  the  wing  surface  has  to  be  generated 
once  (using  the  same  analytic  formulas  used  here  for  the  panel).  The  integrals  are  subsequently  used  for 
wing  box  analysis  and  sensitivities. 

Of  course,  when  the  wing  changes  shape,  panels  change  shape  too.  There  is  linking,  therefore, 
between  overall  wing  planform  design  variables  and  panel  vertex  locations.  This  linking  has  to  be 
accounted  for  in  the  sensitivity  computations  when  overall  wing  /  panel  shapes  are  changing. 


8.7  Eigenvalue  Sensitivities 

The  panel  flutter  problem  was  shown  in  Chapter  2  to  be  a  generalized  eigenvalue  problem  of  the 

form: 


(8.25) 


Let  the  left  eigenvectors  be  {\\f}  from  the  adjoint  eigenvalue  problem: 


(8.26) 


Recall  that  the  matrices 


are  made  up  of  the  mass,  aerodynamic,  stiffness,  and  geometric 


stiffness  matrices  (Eqs.  2.41  and  2.45). 

The  sensitivity  of  the  eigenvalues  with  respect  to  any  design  variable,  DV,  is  obtained  by 
differentiating  Eq.  8.25: 


Premultiplying  by  {\|r}^  gives: 


(8.27) 


-X- 


V 


dDV 


dX  r  -tT 


dDV  dDV 


{Vf  V  {(^} 


(8.28) 
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From  8.26  the  second  term  is  zero  leaving: 


dk 

aov 


’a  If/] 

a  W 

aw 

X  ^  ^  J 

dDV 

'  } 

(8.29) 


The  derivatives  of 


U  and 


are  obtained  using  the  sensitivities  of  all  of  the  system  matrices: 
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ro 

0 

0 

dDV 

ap]_ 

a  K 
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.  dDV 

0 

1  aw 
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dDV 
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a  M 

aov  J 

where 


d  M  _dM 
3  DV  ”  a  DV 


a  C 
aov 


-p.f/, 


aov 


aov  aDv  aov  “aov 


(8.31) 


Therefore,  we  can  calculate  the  eigenvalue  sensitivities  directly  from  the  eigenvectors, 
eigenvalues,  and  the  individual  matrix  sensitivities.  The  eigenvalue  sensitivities  are  useful  to  the  designer 
because  the  real  part  of  the  sensitivity  will  show  the  rate  at  which  a  particular  root  is  going  towards  the 
right-half  s-plane  indicating  how  quickly  that  root  is  nearing  instability. 
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A  more  useful  piece  of  information,  however,  is  the  sensitivity  of  the  flutter  dynamic  pressure 
itself.  If  we  look  at  the  real  part  of  an  eigenvalue,  a,  and  allow  the  design  variable  and  dynamic  pressure 
to  change  then  we  have  a  =  a(qnutter,DV).  At  the  flutter  dynamic  pressure,  the  real  part  of  the  root  is  zero 
by  definition  of  the  flutter  point.  A  variation  of  a  at  flutter  is  then  given  as: 

(8.32) 


which  can  be  rearranged  to  give  the  sensitivity  of  the  flutter  dynamic  pressure  as: 

a  q  flutter  _  8  DV  _  U  DV  J 

aov  da  BX  ^ 

aQflutter  VaQfiu^J 


(8.33) 


The  numerator  is  known  already  from  Eq.  8.29.  The  denominator  is  found  by  replacing  DV  with  qnutter  in 
Eq.  8.29  where 


(8.34) 


The  sensitivity  of  the  flutter  dynamic  pressure  with  respect  to  the  design  variables  tells  the 
designer  where  the  stability  boundary  will  move  if  a  design  variable  is  changed.  The  sensitivity  can  also 
be  used  to  construct  approximations  for  optimization  based  on  approximation  concepts  (Ref.  39).  A  final 
note  on  the  flutter  sensitivity  equation  (Eq.  8.33):  some  special  cases  may  be  encountered  when  multiple 
roots  with  non-distinct  eigenvalues  appear  in  the  stability  analysis.  These  cases  are  beyond  the  scope  of 
this  work,  and  more  information  can  be  found  in  Refs.  40  and  41. 


CHAPTER  9 


VERIFICATION  OF  RESULTS 

9.1  Introduction 

The  numerical  results  of  the  computer  code  implementing  the  techniques  described  in  this  work 
are  compared  with  data  from  the  panel  flutter  literature  to  verify  its  capabilities.  First  the  convergence  of 
the  solution  for  increasing  orders  of  Ritz  polynomials  is  shown  to  establish  the  number  of  terms  necessary 
for  accurate  results.  Then  results  are  presented  for  configurations  of  increasing  complexity  and  compared 
with  available  data.  The  capability  to  accurately  predict  flutter  boundaries  is  shown  for  simply  supported 
skewed  and  trapezoidal  shapes,  yawed  flow,  variable  composite  fiber  orientations,  in-plane  loads,  and 
combinations  of  the  above.  The  analytic  sensitivities  are  compared  with  finite  difference  approximations 
for  validation.  Finally,  the  capability  of  the  code  to  use  analytic  sensitivities  for  direct  and  reciprocal 
approximations  of  flutter  boundaries  is  demonstrated. 

All  of  the  isotropic  test  cases  use  aluminum  material  properties  of  E  =  6.8959  x  10^^  Pa,  pm  = 
2768  kgmass/m^,  and  a  poisson’s  ratio  of  0.3.  All  of  the  composite  test  cases  use  the  material  properties  of 
Ei=137  X  10^  Pa,  E2=9.7  x  10^  Pa,  Gi2=5.5  x  10^  Pa,  pm=1580  kgmass/m^  and  'V)i2=0.3. 


9.2  Convergence  to  Solution 

Results  are  compared  with  Refs.  5  and  23  to  show  the  convergence  rate  of  the  present  analysis. 
The  panels  chosen  from  Ref.  5  are  a  square  and  rectangular  isotropic  aluminum  panel  (Fig.  9.2).  The 
panel  in  Ref.  23  is  a  45®  skewed  panel  made  up  of  a  single  layer  composite  with  a  fiber  angle  of  15®  (Fig. 
9.4).  Figure  9.1  shows  that  satisfactory  convergence  is  obtained  with  either  a  fourth  or  fifth  order 
polynomial  Ritz  series  (Eq.  3.8).  While  the  computational  time  is  significantly  reduced  by  using  a  fourth 
order  polynomial.  Figure  9.1  shows  an  example  of  a  panel  that  requires  a  fifth  order  polynomial.  Most 
panel  configurations  converged  with  fourth  order  series,  but  because  of  the  exceptions,  a  fifth  order 
approximation  is  needed.  Therefore  aU  of  the  results  presented  here  are  based  on  a  fifth  order  Ritz 
polynomial. 


Order  of  Polynomial 


Figure  9.1  Convergence  of  Solution  for  Critical  Dynamic  Pressure  with  no  In-Plane  Loads 


9.3  Isotropic  Rectangular  Panels  in  Yawed  Flow 
Results  are  compared  with  the  analytical  and  FEM  data  of  Ref.  5  to  show  the  effects  of  yawed 
flow  and  rectangular  aspect  ratios.  An  aluminum  plate  was  used  with  the  following  geometry. 

Uoo 


Figure  9.2  Panel  Geometry  for  Ref.  5 
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Reference  5  defines  the  flow  angle  differently  than  the  present  analysis  (Figure  7.1).  The  relation  is  a  = 
-<1).  The  panel  width,  b,  is  held  constant  at  1  .Om  and  panel  thickness,  h,  is  varied  from  2  mm  to  5  mm  so 
that  the  panel  flutters  at  a  meaningful  dynamic  pressure.  The  Mach  number  is  held  constant  at  2.0.  The 
results  are  given  as  a  non-dimensional  critical  dynamic  pressure. 


*  _  ^ 
dVm^  -1 


(9.1) 


Table  9.1  Acrit  for  Rectangular  Isotropic  Panels 
(No  In-Plane  Loads) 


a/b 

Model 

o 

O 

li 

a  =  30‘> 

O 

W-) 

II 

0.5 

Analytical 

385.0 

— 

— 

Sander 

382.0 

213.0 

172.0 

Present 

383.2 

214.9 

175.5 

1.0 

Analytical 

512.6 

— 

— 

Sander 

503.0 

516.0 

523.0 

Present 

511.9 

527.0 

530.7 

1.25 

Analytical 

615.0 

— 

— 

Sander 

612.6 

— 

— 

Present 

614.1 

658 

712 

2.0 

Analytical 

1110 

■  - 

— 

Sander 

1081 

1206 

1388 

Present 

1106 

1201 

1415 

Table  9.1  summarizes  the  results.  It  shows  that  the  present  technique  compares  well  with  both  the 
analytical  and  FEM  results  of  Ref.  5  for  variable  aspect  ratios  and  yaw  angles. 

9.4  Isotropic  Skewed  Panels  Subjected  to  In-Plane  Loads 
We  show  the  effects  of  skewed  panels,  yawed  flow,  and  in-plane  forces  by  comparing  data  with 
the  FEM  results  of  Ref.  21 .  Here  we  use  a  simply  supported  aluminum  panel  with  h  =  2mm  and  =  2.0. 
The  geometry  is  defined  in  Fig.  9.3. 


Figure  9.3  Panel  Geometry  for  Ref.  21 


The  results  are  given  in  the  same  non-dimensional  dynamic  pressure: 


dVm^-1 


(9.2) 


Table  9.2  shows  results  for  a  flow  angle  of  (|)  =  0®  and  Table  9.3  shows  the  same  configurations 
with  (()  =  15"^.  The  present  results  do  not  agree  with  Ref.  21  very  well  for  some  of  the  configurations. 
Recall  from  Chapter  2  that  qcrit  is  found  by  ignoring  damping  and  looking  for  a  coalescence  of 
eigenvalues.  In  most  cases  the  coalescence  frequency  is  between  the  first  and  second  or  second  and  third 
frequencies.  However,  in  the  cases  marked  by  **  our  results  showed  a  coalescence  of  the  14th  and  15th 
roots  before  the  first  and  second  roots.  The  cases  marked  *  had  a  coalescence  of  the  fifth  and  sixth  roots 
before  the  first  and  second  roots.  The  results  of  such  cases  are  shown  as  the  first  and  second  frequency 
coalescence  followed  by  the  higher  pole  coalescence. 


Table  9.2  Acrit  for  Skewed  Panels  with  (j)  =  0° 


us 

Model 

O 

O 

11 

> 

O 

in 

il 

> 

V  =  30® 

0.5 

Ref.  21 

374.1 

303.9 

325.3 

Present 

383.4 

**289/50.7 

242.4 

1.0 

Ref.  21 

518.2 

548.4 

668.2 

Present 

511.9 

518.0 

522.0 

2.0 

Ref.  21 

1147 

1247 

1500 

Present 

1106 

*1118/622 

1138 

Table  9.3  Acrj,  for  Skewed  Panels  with  (j)  =  15 
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L/S 

Model 

o 

O 

II 

> 

O 

II 

> 

1(1  =  30° 

0.5 

Ref.  21 

292.2 

222.1 

243.5 

Present 

274.7 

**222/195 

193.9 

1.0 

Ref.  21 

521.1 

487.0 

579.6 

Present 

515.1 

501.4 

482.6 

2.0 

Ref.  21 

1188 

1169 

1948 

Present 

1134 

*1149/656 

1167 

The  cases  marked  by  *  and  **  in  Tables  9.2  and  9.3  were  re-computed  with  aerodynamic 
damping  for  compmson.  With  the  damping,  the  same  higher  order  poles  converged,  but  there  was 
enough  damping  to  keep  them  stable.  Therefore,  the  model  with  damping  converged  to  flutter  boundaries 
in  better  agreement  with  Ref.  21.  These  additional  test  cases  are  shown  in  Table  9.4.  Reference  30 
discusses  cases  involving  modes  of  nearly  identical  frequencies,  but  weak  aerodynamic  coupling,  that  lead 
to  inaccurate  results  in  the  abscence  of  aerodynamic  damping.  This  may  be  occuring  in  these  cases  and 
should  be  further  investigated. 

Table  9.4  Anmier  and  Acrit  for  Skewed  Panels 


o 

II 

> 

L/S 

Model 

(|)=0° 

<|)=15° 

0.5 

Ref.  21  A<^, 

303.9 

221.1 

Present  Anmt 

306.3 

229.7 

2.0 

Ref.  21  Aorit 

1247 

1169 

Present  Anutt 

1135 

1168 

Finally  we  compare  results  with  in-plane  forces.  Reference  21  uses  a  non-dimensional  in-plane 
force  that  can  be  defined  as: 

N,  1} 


(9.3) 


The  results  are  for  a  square  panel  (L/S=l)  with  (j)  =  0"^. 


Table  9.5  Critical  Dynamic  Pressures  With  ryy  =  rxy  =  0 
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Model 

r„  =  -3 

r«  =  0 

rxx  =  3 

Ref.  21 

275.7 

518.2 

789.0 

Present 

265.1 

512.6 

793.1 

Table  9.6  Critical  Dynamic  Pressures  With  r^x  =  ryy  =  0 


Model 

O 

II 

rxy  —  2 

rxy  =  4 

^xy  “  6 

Ref.  21 

518.2 

487.0 

418.9 

321.5 

Present 

512.6 

473.1 

381.7 

274.3 

Note  the  expected  decrease  in  critical  dynamic  pressure  with  increasing  axial  compression  or  in-plane 
shear. 


9.5  Skewed  Composite  Panels 

To  show  the  effects  of  skew  angle  and  composite  fiber  angle,  results  were  compared  with  Ref  23. 
The  geometry  is  defined  in  Fig,  9.4. 


Figure  9.4  Panel  Geometry  for  Ref.  23 

Table  9.7  compares  the  results  for  the  case  of  an  isotropic  aluminum  panel  with  a=b=l,0m, 
h=3mm  and  Moo=2.0.  The  results  are  in  terms  of  the  nondimensional  critical  dynamic  pressure: 


A. 


(9.4) 


Table  9.7  Acrit  for  Various  Skew  Angles 
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o 

O 

II 

> 

O 

II 

> 

O 

O 

II 

Ref.  23 

511 

522 

664 

Presnt 

513 

532 

606/347* 

With  a  skew  angle  of  30®  the  fifth  and  sixth  poles  coalesced  at  the  lower  dynamic  pressure,  however  this 
did  not  occur  with  aerodynamic  damping. 

Table  9.8  shows  the  critical  dynamic  pressure  for  composite  panels  with  various  fiber  and  skew 
angles.  Table  9.9  shows  the  corresponding  critical  frequencies.  The  results  are  non-dimensionalized  as: 

^crit  = - rrH -  (9-5) 


=co„a^ 


(9.6) 


The  composite  panel  used  had  a=b=l  m  and  h=4  mm. 


Table  9.8  Acrit  for  Skewed  Composite  Panel 


o 

II 

> 

V=30° 

V=45° 

0 

Act  Ref.  23 

Act  Present 

AorRef.23 

Act  Present 

AorRef.23 

Act  Present 

0° 

358 

387 

367 

389 

427 

*429/104 

15° 

248 

246 

299 

301 

321 

319 

O 

O 

170 

163 

202 

195 

257 

254 

45° 

120 

116 

138 

133 

183 

180 

60° 

81 

79.5 

92 

88.5 

160 

128 

75° 

60 

57.5 

87 

73.1 

188 

135 

90° 

69 

57.4 

124 

81.8 

238 

164 

63 


Table  9,9  for  Skewed  Composite  Panels 


\t/=15° 

1 

o 

o 

V=45° 

e 

Qcr  Ref.  23 

Qct  Present 

Ocr  Ref  23 

Qcr  Present 

ftcf  Ref.  23 

ficr  Present 

0“ 

28.4 

30.4 

29.3 

30.9 

34.4 

*34.7/53.3 

15° 

24.3 

TAJ 

30.3 

31.1 

32.6 

33.1 

O 

o 

23.4 

23.5 

26.6 

27.0 

30.7 

31.7 

45° 

22.0 

22.0 

23.0 

23.6 

28.4 

28.6 

60° 

19.3 

19.4 

21.0 

21.0 

30.3 

28.7 

75° 

17.9 

17.9 

22.5 

21.8 

35.8 

33.5 

90° 

19.3 

18.9 

26.6 

24.6 

41.7 

39.0 

The  case  marked  by  *  had  a  coalescence  of  the  third  and  fourth  frequencies  before  the  first  and  second. 
However,  the  results  for  most  of  the  above  configurations  are  very  close  to  Ref.  23.  As  skew  angle  and 
fiber  angle  increase,  correlation  of  the  present  results  and  Ref  23  deteriorates.  Notice  that  the  above 
geometries  are  very  close  to  the  geometries  of  the  skewed  panels  in  Ref  21  but  the  results  compare  much 
more  favorably. 

Data  for  the  flutter  dynamic  pressures  and  frequencies  (aerodynamic  damping  present)  was  also 
obtained  for  these  test  cases  and  is  shown  in  Table  9.10.  The  addition  of  the  aerodynamic  damping 
changed  the  predicted  stability  boundaries  by  only  a  few  percent  and  the  results  are  very  close  to  Ref  23. 
As  before,  with  aerodynamic  damping,  any  higher  frequencies  that  coalesced  remained  stable  and  the 
solution  converged  to  an  appropriate  Anutt. 


Table  9.10  Aflutter  and  fiflutter  for  Skewed  Composite  Panels 


V=15° 

v/=30° 

\|f=45° 

0 

Af  Present 

Qf  Present 

Af  Present 

flf  Present 

Af  Present 

Ilf  Present 

0° 

449 

36.8 

439 

34.1 

447 

35.5 

15° 

268 

26.1 

307 

32.0 

324 

32.2 

o 

O 

166 

23.6 

198 

27.1 

257 

31.8 

45° 

117 

22.0 

134 

23.6 

181 

28.7 

O 

O 

80 

19.4 

89.2 

21.0 

129 

28.7 

75° 

57.8 

17.9 

73.5 

21.8 

135 

33.5 

90° 

57.7 

18.9 

82.3 

24.6 

165 

39.0 

65 


Fig.  9.6  Qcrit  versus  Fiber  Angle  for  Different  Skew  Angles 

The  present  analysis  shows  the  same  effects  of  fiber  angle  on  the  critical  pressure  and  frequency 
as  Ref.  23.  Overall,  except  in  cases  of  high  skew  angle  (45°)  and  high  fiber  angle,  the  results  obtained 
here  are  in  good  correlation  with  the  references. 

9.6  Analytic  Sensitivities 

This  section  shows  the  accuracy  of  the  analytic  sensitivities  and  demonstrates  some  of  the 
capabilities  that  they  give  the  designer.  Figure  9.7  shows  the  accuracy  of  the  analytic  fiber  angle 
sensitivity  compared  with  a  fmite  difference  approximation  for  a  range  of  step  sizes.  The  plot  is  based  on 
the  composite  panel  geometry  shown  in  Fig  9.4  with  v=15°,  0=45°, and  a=b=1.0m.  For  a  step  size 
between  10'^  and  10‘^  the  finite  difference  and  analytic  sensitivities  are  nearly  identical.  A  smaller  step 
size  results  in  an  inaccurate  finite  difference  estimate  because  of  round  off  errors. 


Finite  Difference  Sensitivity/  Analytic  Sensitivi 


The  analytic  sensitivities  above  and  tiie  corresponding  dynamic  pressures  are  used  to  predict  tlie 
effect  of  fiber  angle  on  tlie  stability  boundary  tliat  was  shown  in  secton  9.5.  Figures  9.8  through  9.1 1 
show  tlie  Taylor  series  and  reciprocal  approximations  (Ref  39).  Tliese  approximations,  botli  witli  and 
witliout  aerodynamic  damping,  are  good  for  a  range  of  0  at  least  ±15^, 


Fiber  Angle  (Deg) 


Figure  9.1 1  Anmt  vs.  Fiber  Angle  for  \|f=45'^ 

To  show  the  capabilities  of  the  shape  sensitivities  we  use  Xar  as  an  example.  A  composite  panel 
with  \|/=30°,  0=15'^,  a=b=1.0m,  and  h=3mm  was  used  as  the  referenc  panel.  Xar  was  varied  from  0.5  to 
2.5  m  to  range  from  a  triangular  to  trapezoidal  panel  as  shown  below. 
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Uoo 


\  I 


N 

Figure  9.12  Panel  Geometry  For  Xar  Sensitivity 


tlie  first  order  Taylor  series  approximation  and  the  reciprocal  approximation  give  a  good  estimation  of  the 
affect  of  Xar  on  the  flutter  boundary. 
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Figure  9.13  Anut  vs.  Xar 


Finally,  we  demonstrate  the  capability  of  predicting  the  effect  of  changing  in-plane  loads  on  the 
flutter  boundary.  Recall  that  the  in-plane  load  matrix,  [N],  affects  the  equations  of  motion  linearly 
through  [Kq].  Therefore,  we  can  keep  the  ratio  of  the  Ny’s  the  same  but  allow  for  changing  magnitudes 
by  simply  multiplying  [Kg]  by  a  factor,  t|.  We  used  an  isotropic  panel  defined  by  the  geometry  in  Figure 
9.4  with  \i/=30'^,  a=b=lm,  and  h=3mm.  The  in-plane  loads  were  applied  as  Nx=  -2000,  Ny=  -1000,  and 
Nxy=0.  The  load  factor  was  varied  from  -5.0  to  3.0  to  range  over  tension  and  compression  up  to  the 


^  ^flutt 

dr[ 


-99.99. 


buckling  point.  The  sensitivity  was  computed  at  ti=1.0  to  be 


In-Plane  Load  Factor 


Figure  9.14  Anwt  vs.  In-plane  Load  Factor,  y\ 


The  present  technique  gives  results  in  aggreement  with  the  available  literature  for  a  wide  range 
of  panel  variables.  The  analytic  sensitivities  are  very  accurate  when  compared  with  finite  difference 
approximations,  and  they  have  been  shown  to  successfully  predict  the  effects  of  changing  shape,  fiber 


CHAPTER  10 


CONCLUSION 

The  capability  to  include  configuration  shape  design  variables,  in  any  multidisciplinary  design 
optimization  of  airplanes  in  the  conceptual  or  preliminary  design  stages,  is  essential.  Developments 
in  recent  years  advanced  the  state  of  the  art  in  Design  Oriented  Structural  Analysis  (DOSA)  covering 
approximate  stress,  deformation,  structural  dynamic  and  buckling  analysis  as  well  as  analytical  behavior 
sensitivity  techniques,  where  wing  box  structures,  and  skin  panel  structures,  are  subject  to  shape  as  well  as 
material  and  sizing  variations.  The  present  work  has  focused  on  the  design  oriented  aeroelastic  analysis 
of  optimized  skin  panels  in  supersonic  flow.  Since,  in  typical  optimum  aeroelastic  synthesis  of  wing 
structures,  many  skin  panels  can  be  buckling -critical,  and  since  stressing  a  panel  up  to  a  point  close  to  its 
buckling  load  may  have  a  serious  effect  on  its  aeroelastic  stability  as  a  panel,  it  becomes  important  to 
develop  efficient  analysis  and  sensitivity  capabilities  for  panel  flutter  constraints. 

It  is  shown  in  tliis  work  how  modeling  and  Ritz  formulations  based  on  simple  polynomial 
functions  in  global  coordinates,  lead  to  efficient  evaluation  of  panel  stiffness,  geometric  stiffness,  and 
mass,  as  well  as  aerodynamic  damping  and  aerodynamic  stiffness  matrices.  Using  analytic  formulas  for 
area  integrals  of  polynomial  terms  over  general  trapezoidal  area  shapes,  it  is  shown  that  no  numerical 
integration  is  needed  for  evaluating  panel  structural  or  aerodynamic  matrices.  A  table  of  area  integrals  for 
die  panel,  including  area  integrals  over  the  panel  of  terms  of  the  form  x“  / ,  needs  to  be  evaluated  only 
once  for  a  given  panel  shape.  Then,  structural  and  aerodynamic  matrices,  as  well  as  their  analytic 
sensitivities  with  respect  to  sizing,  material  and  shape  design  variables,  can  be  obtained  by  linear 
combinations  of  elements  of  tliis  table  of  integrals. 

Systematic  evaluation  of  the  resulting  panel  flutter  prediction  capability  was  carried  out, 
comparing  results  from  the  present  work  with  results  from  other  references.  Cases  involving  rectangular 
and  skewed  panels,  isotropic  and  composite  construction,  and  different  types  of  in-plane  loads  were 
covered.  Overall,  the  current  capability  led  to  good  correlation  with  other  prediction  techniques  up  to 
panel  leading  edge  sweep  angles  of  30®.  Large  differences  in  panel  flutter  boundary  predictions  between 
tile  current  technique  and  other  approximate  techniques  including  finite  elements,  were  observed  for 
panel  sweep  angles  of  more  than  30®.  Additional  work  is  needed  in  this  area  to  find  out  whether  the 
discrepancies  are  due  to  limitations  of  the  current  polynomial  Ritz  analysis  or  the  other  approximate  and 
finite  element  techniques  used  for  comparison.  Unfortunately,  only  a  few  papers  on  panel  flutter  address 
the  aeroelastic  stability  of  panels  of  general  trapezoidal  shapes. 

Expressions  for  analytic  sensitivity  of  panel  aeroelastic  poles  and  resulting  flutter  dynamic 
pressure,  have  been  obtained,  and  checked  against  finite  difference  sensitivities.  Excellent  correlation,  and 


a  wide  range  of  step  sizes  for  the  finite  difference  derivatives,  have  been  found.  Used,  in  tuni,  in  direct 
and  reciprocal  Taylor  series  approximations  for  the  flutter  dynamic  pressure,  tlie  flutter  and  sensitivity 
results  have  been  shown  to  lead  to  quite  robust  approximations  over  a  wide  range  of  design  variable 
changes  (wide  move  limits).  The  work  has  also  shown  how  to  integrate  tlie  panel  aeroelastic  analysis  and 
sensitivity  predictions  with  a  wing  box  analysis  and  sensitivity  capability,  where  in-plane  loads 
detennined  by  tlie  wing  box  behavior  serve  as  inputs  to  the  panel  aeroelastic  behavior.  Shape  variations  of 
tlie  wing  and  its  internal  structure  affect  the  panel  both  via  its  in-plane  loads,  and  directly  through  tlie 
effects  of  its  shape  on  its  structural  and  aerodynamic  matrices. 

This  work  was  limited  to  linear  panel  flutter  using  Piston  Theory  aerodynamics,  and  quasi- 
homogeneous  composite  construction.  An  effort  to  extend  the  work  to  cases  of  linear  potential 
aerodynamics,  active  piezoelectric  actuation,  and  non-linear  structural  and  aerodynamic  behavior  is 
currently  underway. 


BIBLIOGRAPHY 


1.  Dowell,  E.H.,  Aeroelasticity  of  Plates  and  Shells,  Noordlioff  International  Publishers,  Leyden,  1975. 

2.  Librescu,  L.,  Elastostatics  and  Kinetics  of  Anisotropic  and  Heterogenous  Shell  Type  Structures, 
Noordhoff  International  Publishers,  Leyden,  1975. 

3.  Dugundji,  J.,  ‘Theoretical  Consideration  of  Panel  Rutter  at  High  Supersonic  Mach  Numbers”,  AIAA 
Journal,  Vol.  4,  No.  7,  July  1966,  pp.  1257-1266. 

4.  Eastep,  F.E.,  and  McIntosh,  S.C.,  “Analysis  of  Nonlinear  Panel  Rutter  and  Response  Under  Random 
Excitation  or  Nonlinear  Aerodynamic  Loading”,  AIAA  Journal,  Vol.  9,  No.  3,  March  1971,  pp.  411-418. 

5.  Sander,  G.,  Bon,  C.,  and  Geradin,  M.,  “Finite  Element  Analysis  of  Supersonic  Panel  Rutter”, 
Intemation  Journal  for  Numerical  Methods  in  Engineering,  7,  pp.  379-394,  1973. 

6.  Yang,  T.Y.,  “Rutter  of  Rat  Finite  Element  Panels  in  Supersonic  Unsteady  Potential  Row”,  AIAA 
Journal,  Vol.  13,  No.  11,  1975,  pp.  1502-1507. 

7.  Bismarck-Nasr,  M.,N.,  “Finite  Element  Analysis  of  Aeroelasticity  of  Plates  and  Shells”,  ASME 
Applied  Mechanics  Reviews,  Vol.  45,  No.  12,  pp.  461-482,  December  1992. 

8.  Ibrahim,  R.A.,  Orono,  P.O.,  and  Madaboosi,  S.R.,  “Stochastic  Rutter  of  a  Panel  Subjected  to  Random 
In  Plane  Forces,  Part  I:  Two  MOde  Interaction”,  AIAA  Journal,  Vol.  28,  No.  4,  April  1990,  pp.  694-702 

9.  Liaw,  D.G.,  “Supersonic  Rutter  of  Laminated  Thin  Plates  with  Thermal  Effects”,  Journal  of  Aircraft, 
Vol.  30,  No.  1,  Jan.-Feb.  1993,  pp.  105-111. 

10.  Shore,  C.P.,  “Rutter  Design  Charts  for  Biaxially  Loaded  Isotropic  Panels”,  Journal  of  Aircraft,  Vol.  7, 
No.  4,  July- August  1970,  pp.  325-329. 

11.  Lin,  K.-J,,  Lu,  P.-J.,  and  Tam,  J.-Q.,  “Rutter  Analysis  of  Composite  Panels  Using  High  Precision 
Finite  Elements”,  Computers  &  Structures,  Vol.  33,  No.  2,  pp.  561-574,  1989. 

12.  Chowdary,  T.V.R.,  Parthan,S.,  and  Sinha,P.K.,  “Finite  Element  Rutter  Analysis  of  Laminated 
Composite  Panels”,  Computers  &  Structures,  Vol.  53,  No.  2,  pp.  245-251, 1994. 

13.  Shiau,  L.-C.,  and  Chang,  J.-T.,  “Transverse  Shear  Effect  on  Rutter  of  Composite  Panels”,  ASCE 
Journal  of  Aerospace  Engineering,  Vol.  5,  No.  4,  October  1992,  pp.  465-479. 

14.  Lee,  I.,  and  Cho,  M.-H.,  “Supersonic  Rutter  Analysis  of  Clamped  Composite  Panels  Using  Shear 
Deformable  Finite  Elements”,  AIAA  Journal,  Vol.  29,  No.  5,  May  1991,  pp.  782-783. 

15.  Hajela,  P.,  and  Glowasky,  R.,  “Application  of  Piezoelectric  Elements  in  Supersonic  Panel  Rutter 
Suppression”,  AIAA  Paper  91-3191,  AIAA  Aircraft  Design  Systems  and  Operations  Meeting,  Baltimore, 
MD,  September  1991. 

16.  Scott,  R.C.,  and  Weishaar,  T.A.,  “Controlling  Panel  Rutter  Using  Adaptive  Materials”,  AIAA  Paper 
91-1007,  April  1991. 


77 

17.  Dowell,  E.H.,  and  Ilgamov,  M„  Studies  in  Nonlinear  Aeroelasticity,  Springer-Verlag,  New  York, 
1988. 

18.  Sipcic,  S.R.,  and  Morino,  L.,  “Dynamic  Behavior  of  Fluttering  Two-Dimensional  Panels  on  an 
Airplane  in  Pull-Up  Maneuver”,  AIAA  Journal,  Vol.  29,  No.  8,  pp,  1304-1312,  August  1991. 

19.  Frampton,  K.D.,  Clark,  R.L.,  and  Dowell,  E.H.,  “State  Space  Modeling  for  Aeroelastic  Panels  Subject 
to  Linearized  Potential  Row  Aerodynamic  Loading”,  AIAA  Paper  95-1294,  Proceedings  of  the  36tli 
AJAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics,  and  Materials  Conference,  New 
Orleans,  LA,  April  10-13,  1995,  Vol.  2,  pp.  1183-1189. 

20.  Durvasula,  S.,  “Rutter  of  Simply  Supported,  Parallelogrammic,  Rat  Panels  in  Supersonic  Row”, 
AIAA  Journal,  Vol.  5,  No.  9,  September  1967,  pp.  1668-1673. 

21.  Kariappa,  Somashekar,  B.R.,  and  Shah,  C.G.,  “Discrete  Element  Approach  to  Rutter  of  Skew  Panels 
witli  IN-Plane  Forces  Under  Yawed  Supersonic  Row”,  AIAA  Journal,  Vol.  8,  No.  11,  November  1970, 
pp.  2017-2022. 

22.  Srinivasan,  R.S.,  and  Babu,  J.C.,  “Free  Vibration  and  Rutter  of  Laminated  Quadrilateral  Plates”, 
Computers  &  Structures,  Vol.  27,  No.  2,  pp.  297-304,  1987. 

23.  Chowdary,  T.V.R.,  Sinha,  P.K.,  and  Parthan,  S.,  “Finite  Element  Rutter  Analysis  of  Composite  Skew 
Panels”,  Computers  &  Structures,  Vol,  58,  No.  3,  pp.  613-620,  1996. 

24.  Haftka,  R.T.,  and  Gurdal,  Z.,  Elements  of  Structural  Optimization,  Kluwer  Academic  Publishers, 
Dordrecht,  Boston,  London,  Third  Edition  1992  (Chapters  6,7) 

25.  Craig,  R.R.,  “Optimization  of  a  Supersonic  Panel  Subjected  to  a  Rutter  Constraint  -  A  Finite  Element 
Solution”,  AIAA  Journal,  March  1979,  pp.  404-405. 

26.  Weisshaar,  T.A.,  “Aeroelastic  Optimization  of  a  Panel  in  High  Mach  Number  Supersonic  Row”, 
Journal  of  Aircraft,  Vol.  9,  1972,  pp.  611-617. 

27.  Weisshaar,  T.A.,  “Panel  Rutter  Optimization  -  A  Refined  Finite  Element  Approach”  International 
Journal  for  Numerical  Methods  in  Engineering,  Vol.  10,  pp.  77-91, 1976. 

28.  Van  Keuren,  G.M.,  and  Eastep,  F.,  “Use  of  Galerkin's  Method  for  Minimum  Weight  Panels  with 
Dynamic  Constraints”,  Journal  of  Spacecraft  and  Rockets,  Vol.  14,  No.  7,  July  1977. 

29.  Librescu,  L.,  and  Beiner,  L.,  “Weight  Mi  nimization  of  Orthotropic  Rat  Panels  Subjected  to  a  Rutter 
Speed  Constraint”,  AIAA  Journal,  Vol.  24,  No.  6,  June  1986,  pp.  991-997. 

30.  Voss,  H.  M.,  and  Dowell,  E.  H.,  “Effect  of  Aerodynamic  Damping  on  Rutter  of  Thin  Panels”,  AIAA 
Journal,  Vol.  2,  No.  1,  January  1964,  ppl  19-120. 

31.  Lottatti,  L,  “The  Role  of  Damping  on  Supersonic  Panel  Flutter”,  AIAA  Journal,  Vol.  23,  No.  10, 
October  1985,  pp.  1640-1642. 

32.  Giles,  G.L.,  “Equivalent  Plate  Analysis  of  Aircraft  Wing  Box  Structures  with  General  Planform 
Geometry”,  Journal  of  Aircraft,  Vol.  23,  No.  11,  November  1986,  pp.  859-864. 


78 

33.  Giles,  G.L„  “Further  Generalization  of  an  Equivalent  Plate  Representation  for  Aircraft  Structural 
Analysis”,  Jounial  of  Aircraft,  Vol.  26,  No.  1,  January  1989,  pp.  67-74. 

34.  Livne,  E„  “Integrated  Multidisciplinary  Optimization  of  Actively  Controlled  Fiber  Composite  Wings”, 
Ph.D.  dissertation,  Mechanical,  Aerospace  and  Nuclear  Engineering  Department,  University  of 
California,  Los  Angeles,  1990. 

35.  Livne,  E.,  “Analytical  Sensitivities  for  Shape  Optimization  in  Equivalent  Plate  Structural  Wing 
Models”,  Journal  of  Aircraft,  Vol.  31,  No.  4,  pp.  961-969,  July-August  1994. 

36.  Livne,  E.,  and  Li,  W-L.,  “Aeroservoelastic  Aspects  of  Wing  /  Control  Surface  Planform  Shape 
Optimization”,  AIAA  Journal,  Vol.  33,  No.  2,  February  1995,  pp,  302-311. 

37.  Livne,  E.,  and  Milosavljevic,  R.,  “Analytic  Sensitivity  and  Approximation  of  Skin  Buckling 
Constraints  in  Wing  Shape  Synthesis”,  Journal  of  Aircraft,  Vol,  32,  No.  5,  September-October  1995,  pp. 
1102-1113. 

38.  Li,  W-L„  and  Livne,  E.,  “Analytic  Sensitivities  and  Approximations  in  Supersonic  and  Subsonic 
Wing  /  Control  Surface  Unsteady  Aerodynamics”,  AIAA  Paper  95-1219,  AIAA/ASME/ASCE/AHS/ASC 
36th  Structures,  Structural  Dynamics  and  Materials  Conference,  New  Orleans,  April  1995. 

39.  Haftka,  R.T.,  and  Gurdal,  Z.,  “Elements  of  Structural  Optimization”,  Third  Edition,  Kluwer  Academic 
Publishers,  Dordrecht,  Boston,  London  1991.  (Chapter  6) 

40.  Pedersen,  P.,  and  Seyranian,  A.,  “Sensitivity  Analysis  for  Problems  of  Dynamic  Stability”, 
International  Journal  of  Solids  and  Structures,  Vol.  19,  No.  4,  pp,  315-335,  1983. 

41.  Seyranian,  A.,  “Sensitivity  Analysis  of  Multiple  Eigenvalues”,  Mechanical  Structures  and  Machines, 
Vol.  21,  No.  2,  pp.  261-284,  1993. 

42.  Golub,  G.,  and  Van  Loan,  C.  F.,  Matrix  Computations,  2nd  edition,  Johns  Hopkins  University  Press, 
1989. 

43.  Schmit,  L.A.,  and  Mehrinfar,  M.,  “Multilevel  Optimum  Design  of  Structures  with  Fiber-Composite 
Stiffened- Panel  Components”,  AIAA  Journal,  Vol.  20,  Number  1,  January  1982,  pp.  138-147. 

44.  Schmit,  L.A.,  and  Ramanathan,  R.K.,  “Multilevel  Approach  to  Minimum  Weight  Design  including 
Buckling  Constraints”,  AIAA  Journal,  Vol.  16,  No.  2,  February  1978,  pp.  97-104. 

45.  Starnes,  J.H.,  Jr,,  and  Haftka,  R.T.,  “Preliminary  Design  of  Composite  Wings  for  Buckling,  Stress  and 
Displacement  Constraints”,  Journal  of  Aircraft,  Vol.  16,  1979,  pp.  564-570. 

46.  Lecina,  G.,  and  Petiau,  C.,  “Advances  in  Optimal  Design  with  Composite  Materials”,  in  Computer 
Aided  Optimal  Design:  Structuyral  and  Mechanical  Systems,  C.A.  Mota  Soares,  editor.  Springer  Verlag, 
Berlin,  1987,  pp.  943-953. 

47.  Ausman,  J.D.,  Hangen,  J.A.,  and  Acevedo,  D.A.,  “Application  of  a  Local  Panel  Buckling  Constraint 
within  Automated  Multidisciplinary  Structural  Analysis  and  Design”,  AIAA  92-1116, 1992  Aerospace 
Design  Conference,  February  3-6, 1992,  Irvine,  CA. 


79 


48.  Ragon,  S.A.,  Gurdal,  Z.,  aiid  Starnes  Jr.,  J.H.,  “Optimization  of  Composite  Box-Beam  Structures 
Including  The  Effects  of  Subcomponent  Interaction”,  AIAA  94-1410,  Proceedings  of  tiie  35th 
AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics  and  Materials  Conference,  Hilton  Kead. 
SC,  April  18-20,  1994,  pp.  818-828. 

49.  Bruhn,  E.F.,  Analysis  and  Design  of  Flight  Vehicle  Structures,  Tri  State  Offset  Company,  1973, 
Chapter  C5. 

50.  Antona,  E.,  and  Romeo,  G.,  “Analytical  and  Experimental  Investigation  on  Advanced  Composite 
Wing  Box  Structures  in  Bending  Including  Effects  of  Initial  Imperfections  and  Crushing  Pressure”,  ICAS- 
86-4.2.3,  15tli  Congress  of  the  International  Council  of  tlie  Aeronautical  Sciences,  London,  England, 
September  1986,  pp.  255-261. 

51.  Jones,  R.M.,  Mechanics  of  Composite  Materials,  McGraw-Hill  Book  Company,  New  York  1975. 

52.  Whitney,  J.  M.,  Structural  Analysis  of  Laminated  Anisotropic  Plates,  Technomic  Publishing 
Company,  Lancaster,  PA,  1987. 

53.  Vinson,  J.  R.,  and  Sierakowski,  R.  L.,  The  Behavior  of  Structures  Composed  of  Composite  Materials, 
Martinus  Nijhoff  Publishers,  1986. 

54.  Tsai,  S.W.,  and  Hahn,  H.T.,  Introduction  to  Composite  Materials,  Technomic  Publishing  Company, 
Lancaster,  PA,  1980. 


